UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

An analysis of lithium-ion battery state-of-health through physical experiments and mathematical modelling Kong, XiangRong (David) 2018

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2018_september_kong_xiangrong.pdf [ 1.86MB ]
JSON: 24-1.0371213.json
JSON-LD: 24-1.0371213-ld.json
RDF/XML (Pretty): 24-1.0371213-rdf.xml
RDF/JSON: 24-1.0371213-rdf.json
Turtle: 24-1.0371213-turtle.txt
N-Triples: 24-1.0371213-rdf-ntriples.txt
Original Record: 24-1.0371213-source.json
Full Text

Full Text

An Analysis of Lithium-ion BatteryState-of-Health through PhysicalExperiments and MathematicalModelingbyXiangRong (David) KongB.Sc., The University of British Columbia, 2016A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2018© XiangRong (David) Kong 2018The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled:   An analysis of lithium-ion battery state-of-health through physical experiments and mathematical modelling       submitted by   XiangRong Kong  in partial fulfillment of the requirements for  the degree of    Master of Science    in    Mathematics     Examining Committee:    Brian Wetton, Mathematics  Co-supervisor   Bhushan Gopaluni, Chemical and Biological Engineering  Co-supervisor Supervisory Committee Member   Arman Bonakdapour, Chemical and Biological Engineering  Additional Examiner     Additional Supervisory Committee Members:      Supervisory Committee Member     Supervisory Committee Member  iiAbstractLithium-ion batteries are ubiquitous in modern society. The high powerand energy density of lithium-ion batteries compared to other forms of elec-trochemical energy storage make them very popular in a wide range of ap-plications, most notably electric vehicles (EVs) and portable devices suchas mobile phones and laptop computers. However, despite the numerousadvantages of lithium-ion batteries over other forms of energy sources, theirperformance and durability still suffer from aging and degradation. The pur-pose of the work presented in this thesis is to investigate how different loadcycle properties affect the cycle life and aging processes of lithium-ion cells.To do so, two approaches are taken: physical experiments and mathematicalmodeling.In the first approach, the cycle life of commercial lithium-ion cells ofLiNiCoAlO2 chemistry was tested using three different current rates to sim-ulate low-, medium-, and high-power consuming applications. The batter-ies are discharged/charged repeatedly under the three conditions, all whiletemperature, voltage, current, and capacity are recorded. Data arising fromthe experiments are then analyzed, with the goal of quantifying batterydegradation based on capacity fade and voltage drop. The results are thenused to build two predictive models to estimate lithium-ion battery state-iiiAbstractof-health (SoH): the decreasing battery V0+ model and the increasing CVcharge capacity model. Furthermore, a simple thermal model fitted fromthe battery temperature profiles is able to predict peak temperature un-der different working conditions, which may be the solution to temperaturesensitive applications such as cellphones.The limitation to physical experiments is that they can be costly andextremely time-consuming. On the other hand, mathematical modeling andsimulation can provide insight, such as the internal states of the battery,that is either impractical or impossible to find using physical experiments.Examples include lithium-ion intercalation and diffusion in electrodes andelectrolytes, various side-reactions, double-layer effects, and lithium concen-tration variations across the electrode layer. Thus, in the second approach,work focuses on implementing the pseudo-two-dimensional (P2D) model, themost widely accepted electrochemical model on lithium-ion batteries. TheP2D model comprises highly-nonlinear, tightly-coupled partial differentialequations that calculate lithium concentration, ionic flux, battery tempera-ture and potential to significant accuracy. The unparalleled prediction abil-ities of the P2D model, however, are shadowed by the low computationalefficiency. Thus, much of this work focuses on reducing model complexity toshorten effective simulation time, allowing for use in applications, such as abattery management system, that have limited computational resources. Inthe end, four model reductions have been identified and successfully imple-mented, with each one achieving a certain standard of accuracy.ivPrefaceThis thesis is original work by the author, XiangRong Kong.A version of Chapter 2 has been published [Kong X.R., Wetton B.,Wilkinson D., Bonakdarpour A., and Gopaluni, B. Advanced Control ofChemical Processes, 2018]. I was the lead investigator, responsible for liter-ary research, experiment design, data analysis, as well as manuscript com-position. Dr. Wilkinson and Dr. Gopaluni from the chemical engineeringdepartment provided the necessary laboratory and equipment to conductexperiments. Dr. Bonakdarpour executed all experiments and manageddata collection. Dr. Wetton was the supervisory author on this project andwas involved throughout the project in concept formation and manuscriptedits.The content of Chapter 3 is based on [18] by the author’s co-supervisor,Dr. Bhushan Gopaluni. However, the author adopted a different approachin solving the original problem, which led to higher efficiency and betterrobustness. Therefore, the content of Chapter 3 is independent of [18] andhence original work by the author.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Lithium-ion Batteries . . . . . . . . . . . . . . . . . . . . . . 21.1.1 Standard lithium-ion battery operation . . . . . . . . 31.1.2 Common lithium-ion battery materials . . . . . . . . 41.2 Battery Degradation . . . . . . . . . . . . . . . . . . . . . . . 41.2.1 Battery use . . . . . . . . . . . . . . . . . . . . . . . . 51.2.2 Rate of charge/discharge . . . . . . . . . . . . . . . . 51.2.3 Temperature . . . . . . . . . . . . . . . . . . . . . . . 5viTable of Contents1.2.4 Depth of charge/discharge . . . . . . . . . . . . . . . 61.3 The Two Approaches . . . . . . . . . . . . . . . . . . . . . . 61.3.1 First Approach: Physical Experiments . . . . . . . . 71.3.2 Second Approach: Mathematical Modelling . . . . . . 71.4 Thesis Structure . . . . . . . . . . . . . . . . . . . . . . . . . 82 State-of-Health Estimation. . . . . . . . . . . . . . . . . . . . . 112.1 Experimental . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.1.1 Galvanostatic cycling of lithium-ion batteries . . . . . 132.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2.1 Voltage vs. Capacity Plots . . . . . . . . . . . . . . . 152.2.2 Temperature Plots . . . . . . . . . . . . . . . . . . . . 162.3 SoH Prediction Model . . . . . . . . . . . . . . . . . . . . . . 182.3.1 SoH prediction model based on decreasing V0+ . . . . 182.3.2 SoH prediction based on increasing CV charging ca-pacity . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.4 Temperature Profile Modelling . . . . . . . . . . . . . . . . . 232.4.1 Temperature profile revisited . . . . . . . . . . . . . . 232.4.2 Temperature fitting . . . . . . . . . . . . . . . . . . . 262.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293 Numerical Modelling of Lithium-Ion Batteries . . . . . . . 303.1 The P2D Model . . . . . . . . . . . . . . . . . . . . . . . . . 303.1.1 Model Overview . . . . . . . . . . . . . . . . . . . . . 303.1.2 Positive and Negative Electrodes . . . . . . . . . . . . 313.1.3 Separator . . . . . . . . . . . . . . . . . . . . . . . . . 36viiTable of Contents3.1.4 Current Collectors . . . . . . . . . . . . . . . . . . . . 373.1.5 Constants and Additional Equations . . . . . . . . . 383.1.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 393.2 Numerical Implementation . . . . . . . . . . . . . . . . . . . 393.2.1 Discretization of Governing Equations . . . . . . . . . 393.2.2 Implementation of Boundary and Interface Conditions 423.2.3 Time Discretization and Newton’s Method . . . . . . 453.2.4 Implementation Results . . . . . . . . . . . . . . . . . 473.2.5 Verification of Implementation Results . . . . . . . . 503.3 Model Reduction . . . . . . . . . . . . . . . . . . . . . . . . . 553.3.1 Two-Parameters Approximation Model . . . . . . . . 573.3.2 Temperature-Reduction Model . . . . . . . . . . . . . 613.3.3 Φs-Reduction Model . . . . . . . . . . . . . . . . . . . 623.3.4 Mixed-Reduction Model . . . . . . . . . . . . . . . . 643.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654 Summary and Future Work . . . . . . . . . . . . . . . . . . . 674.1 First Approach: Physical Experiments . . . . . . . . . . . . . 674.2 Second Approach: Mathematical Modelling . . . . . . . . . . 68Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69AppendicesA Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72viiiList of Tables3.1 First-order convergence in time for Backward Euler method . 533.2 Error and convergence for discretization in x . . . . . . . . . 543.3 Error and convergence for discretization in r . . . . . . . . . . 553.4 Timing comparisons of different models running under 1C,with ∆t = 10, ∆x = 1× 10−6, ∆r = 0.5× 10−6, and h = 1 . . 573.5 Performance of the two-parameters approximation model . . 613.6 Performance of the temperature-reduction model . . . . . . . 623.7 Performance of the Φs-reduction model . . . . . . . . . . . . 633.8 Performance of the mixed-reduction model . . . . . . . . . . . 65ixList of Figures1.1 Schematic of lithium-ion battery being discharged . . . . . . 22.1 Panasonics NCR18650B Lithium-ion Batteries . . . . . . . . 122.2 Discharge/Charge Current and Voltage Profile for 1C Cycling 152.3 Voltage vs. Capacity curves for CCD of a) C/2, b) 1C, c) 2C,and CCC of d) C/2, e) 1C, f) 2C . . . . . . . . . . . . . . . . 162.4 Battery Temperature vs. Cycles . . . . . . . . . . . . . . . . 172.5 V0+ vs. Cycle Number - with Linear Fit, for a) C/2, b) 1C,and c) 2C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.6 Capacity vs. Cycle Number for Constant Voltage Charge -with Linear Fit, for a) C/2, b) 1C, and c) 2C . . . . . . . . . 212.7 Temperature (normalized) vs. Time of CCD and OCV forC/2, 1C, and 2C . . . . . . . . . . . . . . . . . . . . . . . . . 242.8 Average (Vcc − Vd)/2 vs. Capacity for C/2, 1C, and 2C . . . 262.9 Temperature Fitting for CCD and CCC a) C/2 b) 1C c) 2C . 283.1 Visualizing the P2D model . . . . . . . . . . . . . . . . . . . 313.2 Model discretization using the finite difference method . . . . 393.3 One-dimensional finite difference grid structure . . . . . . . . 40xList of Figures3.4 Interpolation technique to recover interface values . . . . . . 433.5 Electrolyte continuity across the cathode and the separator . 443.6 Average temperature of each battery section in 1C dischargewith h = 1W/(m2K) . . . . . . . . . . . . . . . . . . . . . . . 493.7 1C discharge cycle run under different heat exchange coeffi-cients: h = 0.01, h = 1, and h = 100 . . . . . . . . . . . . . . 503.8 Discharge cycle run under 1C, 2C, and 5C . . . . . . . . . . . 513.9 Comparison of the three reduced models with the full model.a) 1C rate comparison. b) 2C rate comparison. c) 5C ratecomparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56xiAcknowledgementsI would like to convey my gratitude to all people who gave me the pos-sibility to complete this thesis. In the first place, I would like to express mysincere gratitude to my supervisor, Professor Brian Wetton, for the contin-uous support of my M.Sc. study and research, for his patience, motivation,and enthusiasm. His logical way of thinking and immense knowledge ofmathematics and physics have been of great value to me. I would also liketo thank my co-supervisor Professor Bhushan Gopaluni, who took muchtime in providing me with valuable ideas and enlightening suggestions, aswell as assisting me with the validation of numerical codes. Moreover, Iwould like to thank Dr. Arman Bonakdarpour for maintaining laboratoryequipment, carrying out battery experiments, and managing experimentaldata, which are crucial to the accomplishment of Chapter 2 of this thesis.Most importantly, I would like to thank my entire family. My immediatefamily to whom this dissertation is dedicated, has been a constant source oflove, concern, support, and strength all these years. I would like to expressmy heart-felt gratitude to my parents who have supported and encouragedme throughout this endeavor. Finally, financial support of the Natural Sci-ences and Engineering Research Council of Canada (NSERC) is gratefullyacknowledged.xiiDedicationTo my wonderful parents, who stood by me and supported all my ideasand dreams. I love you dearly. To my best friends in all aspects of my life,who have been by my side through difficult times. Your constant supportand encouragement continue to help me reaching my goals. I am foreverthankful to have you in my life.xiiiChapter 1IntroductionWith rapidly growing energy consumption, there is a large increase indemand for more efficient and sustainable energy resources. Today, we relyon fossil fuels for most of our energy needs. Combustion of fossil fuels leadsto the emission of greenhouse gases into the atmosphere [7]. Global warmingis a direct consequence of the accumulation of greenhouse gases. Renewablesources such as solar and wind energy are “green” but intermittent energysources. Meanwhile batteries are electrochemical storage devices with whichwe can store energy in the form of chemical potential difference and use itwhenever and wherever needed [6]. During the many decades of research,different chemistries of batteries have been developed, such as nickel cad-mium (NiCd), nickel metal hydride (NiMH), lead acid, lithium ion (Li-ion),and lithium ion polymer (Li-Poly). Among these, lithium-ion batteries pro-vide one of the best tradeoffs in terms of power density, low weight, cellvoltage, and low self-discharge. With these attractive features, lithium-ionbatteries are becoming the uniquitous power sources for portable consumerelectronics such as mobile phones, as well as large battery packs used inhybrid-electric and plug-in hybrid-electric vehicles (HEVs and PHEVs) [11].11.1. Lithium-ion Batteries1.1 Lithium-ion BatteriesA lithium-ion battery comprises four main sections including a positveelectrode (cathode), a negative electrode (anode), an electrolyte, and aporous membrane separator, as shown in Fig. 1.1. The potential differ-ence between the two electrodes drives the lithium ions from one to theother. The current collectors transfer the electrons through the externalcircuit [6].Figure 1.1: Schematic of lithium-ion battery being discharged21.1. Lithium-ion Batteries1.1.1 Standard lithium-ion battery operationThe functioning of a lithium-ion battery is based on the reversible oxida-tion of atomic lithium to lithium ion as given by the following electrochemicalreaction,LiDischargeChargeLi+ + e– (1.1)Before discharge, the positive and negative electrodes, which store lithiumsalts as numerous solid particle grains, are at different chemical potentials.During the discharge process, electrons (e-) are released from the negativeelectrode to the positive electrode through the external circuit. This gen-erates an electrical current to power desired applications. As a result, anelectrical potential difference is developed between the two electrodes. Thelithium-ions (Li+) then detach from the anode solid particles and travel fromthe anode to the cathode through an electrolyte. Finally, the Li-ions inter-calate as atomic lithium in cathode solid particles. The electrolyte, wherelithium-ion migration occurs, is crucial in satisfying the electro-neutralitycondition. The rate at which lithium-ions flow across the interface betweensolid particles and electrolyte is described by the ionic flux j. The transferprocess is reversed when the battery is charging, which is carried out byproviding an external power source to the cell [7].Overall, this redox reaction shown in Eq. 1.1 has a very high electro-chemical potential (Eo = −3.04V compared to the standard hydrogen cell).This property combined with the low molecular weight gives lithium-ionbattery technology an advantage over other alternatives in term of energy31.2. Battery Degradationdensity [11]. The high energy density characteristic has led lithium-ion bat-teries to be popular in a variety of electronic devices. These range from ap-plications with very small energy demand such as implantable cardioverter-defibrillators, to extremely high energy demand such as satellites, and nowhybrid (HEV) and plug-in hybrid (PHEV) electric vehicles [7].1.1.2 Common lithium-ion battery materialsThe commercially available lithium-ion batteries commonly contain graphiteas a negative electrode material, and lithium iron phosphate, lithium cobaltoxide, spinel or their derivatives as a positive electrode material. A typ-ical electrolyte is 1:1 by volume mixture of ethylene-carbonate (EC) anddimethyl-catbonate (DMC) with 1M LiPF6 salt. Poly-olefins have provenstable over time and are universally used in lithium-ion batteries as separa-tors [6].1.2 Battery DegradationDespite the great promise shown by lithium-ion batteries, their perfor-mance suffers from aging and degradation that should be recognized andaccounted for. The degradation of a battery compares to the aging of avehicle fuel tank. A new fuel tank can hold its full capacity and may easilyprovide the engine with all the energy stored at a desired rate. As the tankis refilled, the bottom of the tank may accumulate sediments which slowlydecrease its capacity. These sediments could clog the fuel tank filters whichwould limit the rate at which fuel could be extracted. To make matters41.2. Battery Degradationworse, the fuel tank may develop small holes through which fuel would belost [6]. In the same manner, a battery loses capacity with use throughoutits lifetime, decreasing the total amount of energy it can store and deliver.As the battery degrades, the internal resistance increases resulting in a moredifficult extraction of the stored energy, in addition to the loss of part of thisenergy in the form of heat during charging and discharging. In the followingsubsections some of the causes that may contribute to battery degradationare briefly summarized.1.2.1 Battery useSEI (solid electrolyte interface) protects the negative electrode from ag-gressive solvents in the electrolyte. This film forms during every charge/dischargecycle, however it captures some amount of electrochemically active lithium,increasing the internal resistance and decreasing the capacity.1.2.2 Rate of charge/dischargeHigh rates of charge and discharge are the main contributor for lithiumdeposition, caused when the flow of lithium ions is greater than the onewhich can be intercalated in the electrode. When this occurs, lithium ionsdeposit as metallic lithium which captures active lithium and may even causea short-circuit between the electrodes.1.2.3 TemperatureHigh temperature accelerate unwanted chemical reactions which degradeand reduce the capacity of the battery. With higher temperature, the chem-51.3. The Two Approachesical reaction of reduction on the surface of the carbon electrode will occurat a faster rate, causing gas. With this increased pressure and temperature,the surface film on the electrode can be stretched or damaged. This formscracks through which electrolyte can react with the lithiated carbon par-ticles from the electrode, forming more SEI and therefore capturing moreactive lithium particles.1.2.4 Depth of charge/dischargeA cell which is cycled with high depth-of-discharge (DOD) deterioratesquickly due to a constant stripping and re-depositing of the solid electrodes.A battery will last significantly more cycles with small DODs rather thanwith large DODs. Battery life can be maintained in an acceptable range forplug-in hybrid electric vehicles if very deep cycles, (> 0.6 DOD) are avoided.1.3 The Two ApproachesSignificant research has been focused on trying to understand the ag-ing mechanisms in lithium-ion cells and connect them with measurable andidentifiable features in an effort to improve their utilization and reliability.Understanding why the battery degrades helps to focus efforts where thebattery suffers most. Equally important is the understanding of how muchit degrades [6]. It is crucial for any kind of energy storage to be able toestimate the battery degradation as a function of specific charge/dischargeregimes and conditions. To answer the two questions, there are two popularapproaches: physical experiments and mathematical modeling. This work61.3. The Two Approachestakes on both approaches.1.3.1 First Approach: Physical ExperimentsLithium-ion battery degradation and aging are mainly characterized byirreversible capacity loss (capacity fade) and voltage loss, coupled to changesin kinetic and thermodynamic properties of involved materials and theirinterfaces. In this approach, work centres on the experimental testing oflithium-ion batteries in a laboratory. The batteries are cycled under dif-ferent conditions, all while measurements of temperature, voltage, current,depth of discharge, and cell capacity are performed. Data arising from theexperiments are then scrutinized to form a broad understanding of the cell’saging process, with particular emphasis on quantifying the ability of thecells to store charge as a function of their working life as well as the heatgenerated under different discharge rates. The results are then used to buildtwo predictive models to estimate lithium-ion battery state-of-health (SoH).1.3.2 Second Approach: Mathematical ModellingModelling and simulation can provide insight, such as the internal states(e.g., ionic flux, electrolyte lithium concentration) of the battery that is ei-ther impractical or impossible to find using physical experiments. This ap-proach is particularly important in control applications and battery manage-ment systems (BMS). According to [18], mathematical models for lithium-ion battery dynamics fall within two main categories: Equivalent CircuitModels (ECMs) and Electrochemical Models (EMs). ECMs use only elec-trical components, e.g. inductors, resistors, and capacitors, to represent71.4. Thesis Structurebattery dynamics. While this type of model is structurally simple and com-putationally efficient, cumulative measurement errors, capacity degradationthrough usage life, environmental parameter variation, and device sensitivityto initial conditions heavily affect performance. In contrast, EMs are moreaccurate due to their ability to describe detailed physical phenomena, in-cluding lithium-ion intercalation and diffusion in electrodes and electrolyte,various side-reactions, double-layer effects, and lithium concentration varia-tions. The most widely used EM today is the pseudo-two-dimensional (P2D)model, which is described by a set of tightly coupled and highly nonlinearpartial differential-algebraic equations (PDAEs) [18]. In this approach, wefirst provide a robust implementation of the full P2D model, then discusspotential model reductions to improve simulation efficiency which allows forthe more computationally limiting applications.1.4 Thesis StructureThis thesis is structured into two chapters based on the two approachesdiscussed above.Chapter 2 focuses on the first approach: physical experiments. We firstexplain the type of battery being experimented, the equipment used to per-form discharge/charge cycles, and the tools employed for battery status mea-surements. We then discuss the three treatments (C/2, 1C, and 2C currentrates) used to simulate low-, medium-, and high-power demanding real-lifeapplications as well as the experimental procedures. Next we show batteryaging by providing results on the voltage vs. capacity profiles and also the81.4. Thesis Structuretemperature profiles. Based on the results, two SoH prediction models, thedecreasing V0+ model and the increasing CV charging capacity model, areproposed. Lastly, we take a closer look at the temperature profiles and fitparameters to a simple thermal model, with Newton’s law of cooling as theunderlying thoery. Future work to this chapter involves two projects. First,robustness of the two SoH prediction models can be verified by testing themon different batteries under different conditions. Second, battery aging dueto cycling (discharge/charge) can be distinguished from pure temperaturevariation by dividing batteries into two batches: one is discharged/chargedwith different current rates under room temperature, whereas the other isthermocycled using the temperature profile generated from batch one butnot discharged/charged.Chapter 3 focuses on the second approach: mathematical modeling. Thefirst half of the chapter revolves around the implementation of the full P2Dmodel. We first provide a general overview of the different sections of abattery, then discuss in detail every of the governing P2D equations, as wellas the boundary and interface conditions. Next we define the grid struc-ture and proceed to discretize the entire model using the finite differencemethod. Lastly we implement the model on MATLAB, show results, andanalyze numerical convergence. Although the full P2D model is currentlythe most accurate lithium-ion battery model, it requires an extreme amountof computing resource. However, using a mathematical model for optimiza-tion, parameter estimation, controller design, or life studies requires thesimulation to be run hundreds to thousands of times within a short period,necessitating efficient simulation techinques to be used. Simularly, online91.4. Thesis Structureapplications, such as those found in electric vehicles, have limited compu-tational resources available to do the complex calculations [3]. Therefore,in the second half of the chapter, we discuss potential model reductions to-ward a fast-solver. The four reduced models developed in this chapter, witheach one achieving a different standard of accuracy, become viable optionsin applications that would otherwise be too computationally expensive forthe use of the full P2D model. Future work to this chapter includes codemigration into a compiled language (e.g. C#, Java) for even higher compu-tational efficiency, and also implementation of the reduced P2D models onactual battery controllers. In addition, parallel computing through multi-core CPU, as well as cloud computing on memory-optimized instances canbe experimented to further reduce simulation time.10Chapter 2State-of-Health Estimationof Li-ion BatteriesLithium-ion battery packs are a source of major or supplementary powerfor portable electronics such as cellphones, and laptop computers, mobileapplications such as electric vehicles and electric scooters, and also back-uppower systems of several scales. A key aspect of the technology is theirproprietary Battery Management Systems (BMS) that monitor the batteryto maintain safe operation during charging and use, and allow some per-formance optimization. Such systems have a component that estimates thepack State-of-Charge (SoC), that is, the amount of charge still in the packto deliver power to the application. The simplest SoC indicators rely on aninvariant model of the cell’s performance to yield their output and do nottake into account how a pack is changing over time. However, in reality theperformance of batteries decreases over time and with use, described as achange in the battery’s State-of-Health (SoH). Specifically, the battery SoHis defined by a lower battery capacity and an increased internal resistance.In recent years, a lot of attention has been focused on the diagnosis oflithium-ion battery SoH. Nearly all literature regard discharge/charge cur-112.1. Experimentalrent rate and operation temperature as top factors affecting SoH. However,a concrete mathematical model which can connect these two factors to bat-tery’s SoH is still not available.This chapter introduces experiments that involve discharging and charg-ing lithium-ion batteries over many cycles, while their voltage, capacity,current, and temperature profiles are recorded. Based on the data analysis,two predictive lithium-ion battery SoH models are presented.2.1 ExperimentalCylindrical 18650 lithium-ion rechargeable cells (Panasonic NCR18650B,Fig. 2.1) of Lithium Nickel Cobalt Aluminum Oxide (LiNiCoAlO2) chem-istry are tested in this work. The nominal cell voltage and capacity are 3.6Vand 3.2Ah, respectively. The manufacturer recommended charge/dischargevoltage boundaries are between 2.5 and 4.2V.Figure 2.1: Panasonics NCR18650B Lithium-ion Batteries122.1. Experimental2.1.1 Galvanostatic cycling of lithium-ion batteriesGalvanostatic cycling refers to the discharge and charge of batteries with“static”, or constant, current rates. It is performed under room tempera-ture within the manufacturer specified voltage range of 2.5-4.2V. Each cycleconsists of six stages,1. Constant current discharge (CCD) at a chosen C-rate until 2.5V2. Open circuit voltage (OCV) for 30 minutes3. Constant current charge (CCC) at 1C until voltage reaches 4.2V4. OCV for 30 minutes5. Constant voltage charge (CVC) at 4.2V for 4 hours6. OCV for 30 minutesHere, the C-rate is the current rate required to discharge the battery within1 hour to the manufacturer specified voltage cutoff. In our case, the voltagecutoff is 2.5V and 1C=3.2A. Meanwhile, the open circuit voltage (OCV) insteps (2), (4) and (6) refers to disconnecting any external circuit so thatvoltage and electrolyte inside the battery equilibrate. In this work, thechosen C-rates for step (1) are C/2, 1C, and 2C, in order to simulate low-,medium-, and high-power demanding applications. Note that all cycles wereperformed under room temperature, ∼ 22o Celsius.A 10A, VMP3 multi-channel potentiostat from BioLogic Science Instru-ments was used to perform the discharge/charge cycles, with 30 minutes ofOCV between each discharge/charge step, all while voltage, current, and132.2. Resultscapacity are recorded in 30 second intervals. Temperature of the cell duringcycling was recorded by self-adhesive, silicon based cement, K-type thermo-couple temperature sensors (SA1K-72, Omega Engineering Inc.) with betterthan 0.3s response time. Two thermocouples were used: one attached to themain body of the cell and one attached to the anode tip. The cover of thecell was carefully removed where the self-adhesive thermocouple was at-tached. Temperature data were recorded by a high-speed 8-channel TCICthermocouple interface card (TCIC-USB-ENC, Omega Engineering Inc.).Fig. 2.2 shows a typical 1C current and voltage profile of the battery forone discharge/charge cycle. During constant current discharge, the currentis denoted as a negative value (-3.2A) and the voltage decreases nonlinearlyuntil it reaches the cutoff voltage of 2.5V. During constant current charge,the current is held constant at 3.2A with a rising voltage until the voltagereaches the maximum voltage of 4.2 V, at which point constant voltagecharge begins. The profiles for C/2 and 2C follow the same pattern exceptthe current rates are 1.6A and 6.4V, respectively.2.2 ResultsRecall that each discharge cycle completes in approximately i) 2hrs forC/2, ii) 1hr for 1C, and iii) 0.5hr for 2C, excluding the constant currentand constant voltage charge processes and the OCVs in between. Therefore,due to the limitation of time and equipment, 200 cycles of 2C, 100 cycles of1C, and 25 cycles of C/2 are performed.142.2. ResultsFigure 2.2: Discharge/Charge Current and Voltage Profile for 1C Cycling2.2.1 Voltage vs. Capacity PlotsVoltage vs. capacity curves give crude estimation of battery degradation.As a battery is cycled, its maximum capacity decreases. This is indicatedby a leftward shift of the curves.Fig. 2.3 shows voltage vs. capacity curves. From top to bottom, eachrow represents C/2, 1C, and 2C current rate, respectively. The left columnrepresents constant current discharge whereas the right column representsconstant current charge. Notice that both discharge and charge curves of allthree current rates illustrate a similar pattern: they move to the left overtime, indicating that the capacity, or the amount of charge the battery isable to hold, decreases as it is cycled. Going down the rows, we observe thatas C-rate increases, the capacity drop is on average larger every cycle, and152.2. Resultsalso the starting voltage of the discharge is lower. Both of these illustratethat higher current rates cause faster battery aging.Figure 2.3: Voltage vs. Capacity curves for CCD of a) C/2, b) 1C, c) 2C,and CCC of d) C/2, e) 1C, f) 2C2.2.2 Temperature PlotsFig. 2.4 shows the battery temperature vs. time of each cycle for C/2,1C, and 2C. As labeled on the plots, the first, second, and third peak eachrepresents the constant current discharge, constant current charge, and con-stant voltage charge, respectively, while the three temperature drops rep-162.2. Resultsresent the OCVs. Notice that there are relatively large room temperaturefluctuations during the cycling process. Going from top to bottom, we ob-serve that higher current rates cause faster and larger rise in battery temper-ature. However, there is no clear association between battery temperatureand cycle number.Figure 2.4: Battery Temperature vs. Cycles172.3. SoH Prediction Model2.3 SoH Prediction ModelPresent models of lithium-ion battery SoH prediction are based on largeamounts of battery operation data. However, even the most advanced bat-tery management system (BMS) nowadays has very limited data storage.Based on the assumption that voltage vs. capacity curves follow a simi-lar pattern for varying discharge rates, the models presented in this sectionattempt to predict SoH through only a short history of battery operation,which makes BMS implementation highly feasible. Specifically, given merelythe voltage and charge drawn, the two models predict the number of dis-charge/charge cycles the battery has performed through, a natural indicatorof battery SoH.2.3.1 SoH prediction model based on decreasing V0+Recall that as the battery is cycled, the amount of charge it is able tostore decreases, resulting in a loss of power output. Let γ denote the cyclenumber a battery has been discharged/charged through. Let V0+ be thevoltage at the beginning of discharge (see Fig. 2.2), with it a function ofcycle number,V0+ = V0+(γ) (2.1)We expect V0+ to decrease as battery ages. This property alone can be anindicator of battery SoH. Fig. 2.5 shows the plot of V0+ vs. γ for the threeC-rates.Notice that the vertical axis for the three plots are intentionally left un-normalized for the best visual inspection. All three plots show that V0+ and182.3. SoH Prediction ModelFigure 2.5: V0+ vs. Cycle Number - with Linear Fit, for a) C/2, b) 1C, andc) 2Cγ is linearly related. Performing linear regression, we obtain the function,V0+(γ) = β0 + β1γ (2.2)The linear fit superimposed on the data points is shown on the same figure.The values of β0 and β1 for the three C-rates are listed in the table below.192.3. SoH Prediction ModelC/2 1C 2Cβ0 4.0062 3.7831 3.4928β1 −0.00056 −0.00068 −0.0013We can interpret β0 as the V0+, in Volts, of a new battery which has notgone through any discharge/charge cycles, and β1 as the rate of decrease inV0+, in Volts per cycle. Since the batteries were new and never used prior tothe experiments, we observe that β0 decreases with increased current rate asexpected, while β1 drops. This make sense as higher current rate results infaster battery aging, indicated by a faster rate of decrease in V0+ per cycle.With β0 and β1 found for the three C-rates, to predict lithium-ion bat-tery SoH using this model, simply determine the C-rate of the application,and then evaluate cycle number γ using Eq. 2.2, with the battery’s V0+measured. A small γ implies a relatively healthy battery whereas a large γimplies the opposite. The limitations to this model are that it assumes astrictly constant current rate of the battery’s application, and also it requiresthe battery to be fully charged to access V0+.2.3.2 SoH prediction based on increasing CV chargingcapacityThe previous SoH prediction model assumes the current rate of dis-charge/charge process is constant, which is usually not the case in real ap-plications. For example, the driver of an electric vehicle may accelerate fora few seconds, drawing a large current, followed by an immediate deceler-ation, drawing a small current. However, the constant voltage charge after202.3. SoH Prediction Modelthe constant current charge is relatively much more consistent, since it doesnot depend on the previous discharge/charge current rates. Hence, measur-ing the change in capacity during constant voltage charge as a function ofcycle number is another way to assess the SoH of the battery.Figure 2.6: Capacity vs. Cycle Number for Constant Voltage Charge - withLinear Fit, for a) C/2, b) 1C, and c) 2C212.3. SoH Prediction ModelFig. 2.6 shows the capacity vs. cycle number of constant voltage chargefor the three C-rates. A linear fit of the data points is superimposed on thesame figure, with the equation,capacity(γ) = δ0 + δ1γ (2.3)We observe that, independent of the C-rate, the amount of charge neededto complete the constant voltage charge process is positively related to thecycle number. This is expected since as battery ages, its internal resistanceincreases due to SEI formation, so charging is more difficult and takes longer.Going down the rows, as C-rate increases, the slopes of capacity vs. cyclenumber plots become steeper, implying faster battery aging. The values ofδ0 and δ1 for the three C-rates are listed in the table below.C/2 1C 2Cδ0 1.1652 1.1691 1.1723δ1 0.0002 0.0014 0.0016To predict lithium-ion battery SoH using this model, simply measurethe battery capacity during constant voltage charge and then evaluate cyclenumber γ using Eq. 2.3. Compared to the previous model, this one has theadvantage that it does not require the application of the battery to drawa strictly constant current. The downside is that measuring V0+ for theprevious model takes a single data point whereas measuring capacity takesmany, depending on the granularity of measurement. Therefore this modelrequires the BMS to have a more advanced data storage. Furthermore,notice that the parameters δ0 and δ1 for 1C and 2C are close in magnitude,222.4. Temperature Profile Modellingimplying that SoH prediction may have a lower accuracy for high currentrates.2.4 Temperature Profile ModellingWe now analyze the temperature profile recorded and fit parameters toa simple thermal model.2.4.1 Temperature profile revisitedThe column on the left on Fig. 2.7 shows the close-up of the constantcurrent discharge and OCV of the temperature profile of each cycle for thethree C-rates. Here the rising temperature is caused by constant currentdischarge, while the temperature drop to ambient represents the OCV. Thelinear trend in temperature has been removed to reduce the effect of roomtemperature variations.The temperature profiles for all cycles are averaged and plotted in theright column of Fig. 2.7. The decrease of the temperature is in an exponen-tial decay form, suggesting Newton’s law of cooling, which can be modelledby,∂T∂t= −k(T − Troom) (2.4)where k is a constant to be fitted, and Troom is the room temperature.The rising temperature for every C-rate, on the other hand, begins toshow a plateau after midway (expected if the heat generation was constant intime), but experiences an inflection point, as labeled on the plots, and thenincreases sharply until constant current charge is over. This phenomenon232.4. Temperature Profile Modellingis obvious for C/2 and 1C, and can also be observed for 2C under a morecareful scrutiny. This suggests that there is an increasing volumetric heatgeneration term, g(t), during discharge,∂T∂t= −k(T − Troom) + cg(t) (2.5)where k is the same constant as in Eq. 2.4, and c is another constant to befitted. The inverse of c is the average thermal capacity of the cell.Figure 2.7: Temperature (normalized) vs. Time of CCD and OCV for C/2,1C, and 2C242.4. Temperature Profile ModellingAccording to [18], this volumetric heat generation term has the form,g(t) =I(IR+ η)pir2H(2.6)where I is the current, R is the effective resistance, η is the overpotential,and r and H are radius and height of the cell, respectively. The denominatorrepresents the volume of the cylindrical cell.Decoupling the resistance from the overpotential is challenging from thisset of data. We consider the average of the voltage difference between chargeand discharge for C/2, 1C, and 2C in Fig. 2.8. The difference, divided by2, is relatively constant in discharge capacity initially. Thus we consider theinitial (Vcc − Vd)/2I as the effective “resistance” R. Note that in Fig. 2.8,the discharge capacity for 1C begins at 0.4: as the cell charges, constantvoltage charge begins when the cell voltage reaches 4.2V and the remainingcharge is added in a way that is not directly comparable to discharge. Thesame pattern applies to the other two C-rates.Although there is an overpotential for both charge and discharge, weascribe the overpotential η in Eq. 2.5 to predominantly the discharge, thus,IR+ η =(Vcc − Vd)∗/2 θ < θ∗Vcc − Vd − (Vcc − Vd)∗/2 θ > θ∗where θ∗ is the capacity at which constant current charge changes to constantvoltage charge (θ∗ ≈ 1.2Ah for 1C as shown in Fig. 2.8), and (Vcc − Vd)∗/2is the value of the voltage difference at this capacity (≈ 0.40V in this case).252.4. Temperature Profile ModellingFigure 2.8: Average (Vcc − Vd)/2 vs. Capacity for C/2, 1C, and 2C2.4.2 Temperature fittingThe analytical solution to Eq. 2.4 is,T (t) = Ti + (Tm − Ti)ek(tm−t) (2.7)and the solution to Eq. 2.5 is,T (t) = Ti + c∫ t0g(s)ek(s−t)ds (2.8)where ti and Ti are respectively starting time and starting temperature ofthe temperature rise, and tm and Tm are respectively starting time and262.4. Temperature Profile Modellingstarting temperature of the temperature drop.By fitting numerical solutions of Eq. 2.7 and 2.8 to experimental dataand then minimizing residuals, the values of k and the inverse of c for eachC-rate are obtained as,C/2 1C 2Ck 0.0180 0.0204 0.01941/c 760.17 731.31 735.04The quality of the fit for the three C-rates is shown in Fig. 2.9. Bydimensional analysis, k has units of 1s and c has units ofJm3K. From theresults, k is simply the constant in Newton’s law of cooling while c is theinverse of average thermal capacity [18].From Eq. 2.7, we notice that the term (tm − t) in the exponent isalways negative. Therefore, the value of k determines the rate in whichbattery temperature equalizes with room temperature: a larger k impliesa faster rate and a smaller k implies the opposite. Therefore, k representstemperature insulation with the outside, an intrinsic property of the batterymaterial. Since the batteries tested in the experiments are of the exact sametype and thus of the exact same material, we observe from the table thatthe k are similar in value and independent of the current rate. On the otherhand, the 1/c term represents the amount of energy in joules needed toincrease the battery temperature by 1o. This is also an intrinsic propertyof the battery and is independent of the current rate applied. Hence, weobserve that the values of 1/c are similar across C/2, 1C, and 2C, as shownin the table.272.4. Temperature Profile ModellingFigure 2.9: Temperature Fitting for CCD and CCC a) C/2 b) 1C c) 2CAn useful application the simple thermal model in this section providesis a predictive algorithm of battery temperature. With voltage, current, andtemperature profiles measured, this thermal model can provide informationon intrinsic properties of the battery as represented by k and 1/c. Thisinformation can then be used to answer questions on the peak temperaturecertain applications will reach under given environmental temperature andcurrent rate. The answers to these questions ultimately concerns the mostcrucial problem of battery usage: is the battery safe within the predictedtemperature range.282.5. Conclusion2.5 ConclusionIn this study, Panasonic NCR18650B lithium-ion batteries were cycledunder room temperature with current rates of C/2, 1C, and 2C while theircurrent, voltage, capacity, and temperature profiles were recorded. Bat-tery degradation and aging is clear from the drop of V0+ and the leftwardshift of voltage vs. capacity curves. Two models were proposed to predictlithium-ion battery state-of-health: 1) the decreasing battery V0+ modeland 2) the increasing CV charge capacity model. Each has its advantageand limitations. Additionally, we derived a fitted thermal model that canbe used to predict cell temperature under different current rates, which isuseful if the container of the battery is temperature sensitive, for example,a mobile phone. Future work that may yield interesting results includesdistinguishing battery aging due to cycling (discharge/charge) from purelytemperature variation. This can be done by setting up two treatments: thefirst treatment involves discharging/charging batteries with different cur-rent rates under room temperature, while the second treatment involves nodischarge/charge but thermocycles batteries using the temperature profilegenerated from the first treatment.29Chapter 3Numerical Modelling ofLithium-Ion Batteries3.1 The P2D Model3.1.1 Model OverviewThe pseudo-two-dimensional (P2D) model consists of coupled nonlinearPDAEs for the conservation of mass and charge in the three sections of thebattery - cathode, separator, and anode - denoted respectively by the in-dexes p, s, and n. The positive and negative current collectors are denotedby a and z. The index i ∈ S is used to refer to a particular section of the bat-tery, where S := {a, p, s, n, z}. The governing equations of the P2D modelare taken from [18]. The objectives of this work are to (i) provide a robustimplementation, (ii) identify model reductions, and (iii) make progress to-wards a fast-solver. Fig. 3.1 depicts the five domains inside of the batterycell as well as the virtualization of solid particles inside the two electrodes.303.1. The P2D ModelFigure 3.1: Visualizing the P2D model3.1.2 Positive and Negative ElectrodesSolid-particle concentrationDiffusion inside solid spherical particles with radius Rp is described byFick’s law,∂cs(r, t)∂t=1r2∂∂r[r2Dsp∂cs(r, t)∂r](3.1)313.1. The P2D Modelwith boundary conditions,∂cs(r, t)∂r∣∣∣r=0= 0,∂cs(r, t)∂r∣∣∣r=Rp= −j(x, t)Dseff(3.2)where r is the radial direction, or the pseudo-second-dimension, along whichthe ions intercalate within the active particles. Here j represents the ionicflux across the solid particles and the electrolyte.Solid-particle potentialSolid-particle potential in the two electrodes, Φs(x, t) ∈ R, is describedby the equation,∂∂x[σeff,i∂Φs(x, t)∂x]= aiFj(x, t) (3.3)Due to physical constraints, it is necessary to impose zero-flux boundaryconditions for Φs at the interface between electrodes and the separator, aswell as the enforcement of Ohm’s law at the cathode and anode ends,σeff,i∂Φs(x, t)∂x∣∣∣x=xˆ0,xˆn= −Iapp(t) (3.4)σeff,i∂Φs(x, t)∂x∣∣∣x=xˆp,xˆs= 0 (3.5)Here Iapp(t) is the applied current density given as an operating condition.323.1. The P2D ModelElectrolyte concentrationIn the positive and negative electrodes, the electrolyte concentrationce(x, t) ∈ R+ is described by the equation,i∂ce(x, t)∂t=∂∂x[Deff,i∂ce(x, t)∂x]+ ai(1− t+)j(x, t) (3.6)where t ∈ R+ is the time and x ∈ R is the spatial direction through electrodesand separator along which the ions are transported. The first term on theright represents diffusion of the electrolyte while the second term representsionic flux from the solid particles.At the cathode and anode ends, we impose zero-flux boundary condi-tions,∂ce∂x∣∣∣x=x0,xn= 0 (3.7)Meanwhile at the two electrode-separator interfaces, we enforced the conti-nuity of electrolyte concentration,ce(x, t)∣∣∣x=x−p= ce(x, t)∣∣∣x=x+p(3.8)ce(x, t)∣∣∣x=x−s= ce(x, t)∣∣∣x=x+s(3.9)Similarly, continuity of fluxes is also enforced. Due to changes in materialproperties along the length of the battery, the values of different coefficients(e.g., Deff,i, κeff,i, λi) need to be evaluated at the interface between twodifferent materials. For the flux of electrolyte at the two electrode-separator333.1. The P2D Modelinterfaces, we have−Deff,p∂ce(x, t)∂x∣∣∣x=xˆ−p= −Deff,s∂ce(x, t)∂x∣∣∣x=xˆ+p(3.10)−Deff,s∂ce(x, t)∂x∣∣∣x=xˆ−s= −Deff,n∂ce(x, t)∂x∣∣∣x=xˆ+n(3.11)Electrolyte potentialElectrolyte potential in the two electrodes, Φe(x, t), is described by theequation,aiFj(x, t) = − ∂∂x[κeff,i∂Φe(x, t)∂x]+∂∂x[κeff,iΥT (x, t)∂ ln ce(x, t)∂x](3.12)Given that only potential differences are measurable, without loss of gener-ality, Φe can be set to zero at the end of the anode. On the cathode side,zero-flux conditions are imposed,∂Φe∂x∣∣∣x=x0= 0, Φe∣∣∣x=xn= 0 (3.13)At the two electrode-separator interfaces, similar to the electrolyte concen-tration, continuity of potential,Φe(x, t)∣∣∣x=x−p= Φe(x, t)∣∣∣x=x+p(3.14)Φe(x, t)∣∣∣x=x−s= Φe(x, t)∣∣∣x=x+s(3.15)343.1. The P2D Modelas well as continuity of fluxes,−κeff,p∂Φe(x, t)∂x∣∣∣x=x−p= −κeff,s∂Φe(x, t)∂x∣∣∣x=x+p(3.16)−κeff,s∂Φe(x, t)∂x∣∣∣x=x−s= −κeff,n∂Φe(x, t)∂x∣∣∣x=x+s(3.17)are enforced.TemperatureTemperature variations are also included with the set of equations de-scribing the system. The thermal equations include different source terms,which are the ohmic, reversible, and reaction generation rates Qohm, Qrev,and Qrxn, respectively.,ρiCp,i∂T (x, t)∂t=∂∂x[λi∂T (x, t)∂x]+Qohm +Qrxn +Qrev (3.18)The ohmic generation rate takes into account heat generated as a conse-quence of the motion of lithium-ions in the solid/liquid phase. The re-action generation rate accounts for heat generated due to ionic flux andover-potentials, and the reversible generation rate takes into account theheat rise due to the entropy change in the electrodes’ structure [18].At all section interfaces, boundary conditions include both continuityof solution and continuity of flux. For example, at the cathode-separator353.1. The P2D Modelinterface,−λp∂T (x, t)∂x∣∣∣x=x−p=− λs∂T (x, t)∂x∣∣∣x=x+p(3.19)T (x, t)∣∣∣x=x−p=T (x, t)∣∣∣x=x+p(3.20)Boundary conditions at other interfaces are of similar form.Ionic FluxIntertwining temperature, electrolyte concentration, electrolyte poten-tial, solid-particle concentration, and solid-particle potential is j(x, t). j(x, t)is the flux of lithium ions across the surface of the solid-particles into theelectrolyte at position x and time t, and is given by Butler-Volmer kinetics,j(x, t) =2keff,i√ce(x, t)(cmaxs,i − c∗s(x, t))c∗s(x, t) sinh[ 0.5FRT (x, t)ηi(x, t)](3.21)whereηi(x, t) =Φs(x, t)− Φe(x, t)−Ui (3.22)represents the overpotential. Note that i = {p, n}, indicating the ionic fluxis present in only the positive and negative electrodes but not the separator.3.1.3 SeparatorSince the separator is absent of any solid particles, the dynamics inthe separator is simplified as equations of solid-particle concentration and363.1. The P2D Modelpotential, cs(r, t) and Φs(x, t), as well as the ionic flux, j(x, t), are eliminated.Electrolyte concentrationThe ce equation of the separator, in contrast to that of the electrodes,consists of purely diffusion and no ionic flux,i∂ce(x, t)∂t=∂∂x[Deff,i∂ce(x, t)∂x](3.23)Similarly, the electrolyte potential is also independent of the ionic flux,0 = − ∂∂x[κeff,i∂Φe(x, t)∂x]+∂∂x[κeff,iΥT (x, t)∂ ln ce(x, t)∂x](3.24)3.1.4 Current CollectorsThe two current collectors span the two ends of the battery. Absent ofboth electrolyte and solid particles, temperature rise in the current collectorsis caused solely by the applied current density,ρiCp,i∂T (x, t)∂t=∂∂x[λi∂T (x, t)∂x]+I2app(t)σeff,i, i ∈ a, z (3.25)and Newton’s law of cooling with the outside,− λa∂T (x, t)∂x∣∣∣x=0= h(T ref − T (x, t)) (3.26)− λz ∂T (x, t)∂x∣∣∣x=L= h(T (x, t)− T ref) (3.27)373.1. The P2D ModelThe heat exchange coefficient h is proportional to the reciprocal of temper-ature insulation: a low h indicates high insulation and faster increase ofbattery temperature, and the opposite for high h.3.1.5 Constants and Additional EquationsAll experimentally measured parameters and additional equations usedin the implementation are reported in Table I and II in the Appendix. Theopen circuit voltage (OCV) is denoted by U while the entropic variation ofthe OCV is denoted by ∂U∂T . Since the cathode, anode, and the separator arecomposed of different materials, for a given section i, different electrolytediffusion coefficients Di, solid-phase diffusion coefficients Dsi , electrolyte con-ductivities κi, porosities i, thermal capacities Cp,i, thermal conductivitiesλi, densities ρi, solid-phase conductivities σi, particle surface area to vol-umes ai, maximum solid phase concentrations cmaxs,i , overpotentials ηi, andparticle radiuses Rp,i, can be defined. The terms R and F are the universalgas constant and the Faraday constant, repsectively, with t+ representingthe transference number and Tref the environment temperature.Within the battery, continuous interface conditions are imposed acrossthe different materials. In order to get a more detailed description of theconductivity (κeff,i) and diffusion phenomena (Deff,i) inside the electrolyte,all the related coefficients are determined as a function of ce and T . In orderto take into account the properties of different materials used in the battery,effective diffusion and conductivity coefficients are evaluated according tothe Bruggeman’s theory, with “eff” suffixes representing effective values ofeach coefficients.383.2. Numerical Implementation3.1.6 SummaryIn this section, all governing equations, interface and boundary condi-tions of the P2D model are explained. Mathematically, the system is mixedof second-order boundary value and parabolic type. This is a typical struc-ture of models coming from electrochemical systems [9]. We are now readyto proceed to the model implementation.3.2 Numerical Implementation3.2.1 Discretization of Governing EquationsRecall that the battery is composed of five sections: positive currentcollector (a), cathode (p), separator (s), anode (n), and negative currentcollector (z). The cathode and the anode each further contains solid spher-ical particles with radius Rp, resulting in the pseudo-second dimension r.The overall picture of the model is depicted in Fig. 3.2.Figure 3.2: Model discretization using the finite difference methodDimension x and pseudo-second-dimension r are both discretized on astaggered grid using the finite difference method. The grid structure inthe x-direction is defined by subdividing the spatial domain x ∈ R into393.2. Numerical ImplementationNa+Np+Ns+Nn+Nz non-overlapping segments with geometrically centerednodes (as depicted in Fig. 3.3). Every segment is associated with a centre xnand spans the interval [xn− 12, xn+ 12]. The unknown variable at xn is denotedby Ωn.Figure 3.3: One-dimensional finite difference grid structureTo facilitate the treatment of boundary and interface conditions, theends of each segment are aligned with the domain boundaries and internalinterfaces. The number of segments in each section, Ni for i ∈ {a, p, s, n, z},is chosen so that the width of every segment is uniform across all five sectionsand is defined as∆x =∑i li∑iNi(3.28)where li represents the length of a particular section of the battery and islisted in Table IV.At each xn, the pseudo-second-dimension r is discretized using the sameapproach except it is only present in the cathode and anode and has adifferent segment width ∆r.Once the discretization grid is structured, the governing equations arediscretized with finite difference. The central difference scheme is used forboth first and second derivatives. A few key discretizations are shown in this403.2. Numerical Implementationsections. All the interface conditions used to enforce continuity between ad-jacent materials are discussed in Implementation of Boundary and InterfaceConditions section.Solid-particle concentrationThe solid-particle concentration equation is discretized as follows∂cs,n(r, t)∂t=Dspr2mr2m+1/2(cs,n,m+1 − cs,n,m)− r2m−1/2(cs,n,m − cs,n,m−1)∆r2(3.29)where rm is the coordinate of dimension r measured from the center of theparticle. The solid-particle surface concentration, c∗s, which is needed in theionic flux equation, can be obtained using the ghost point technique,c∗s,n =cs,n,M + cs,n,M+12(3.30)where suffixes s,n,M and s,n,M+1 represent the last and the ghost point ofthe solid particle at a particular x.Solid potentialThe solid potential equation of the electrodes is discretized as follows,σeff,iΦs,n−1 − 2Φs,n + Φs,n+1∆x= aiFjn (3.31)In the separator, since there is no solid particles, this equation does notapply.413.2. Numerical ImplementationTemperatureFor the temperature equation in the electrodes, the reversible and reac-tive heat sources can be discretized asQrxn,n =Faijnηi,n (3.32)Qrev,n =FaijnTn∂Ui,n∂T(3.33)whereas the derivatives present in the ohmic source are numerically approx-imated as∂Φs(x, t)∂x|xn ≈Φs,n+1(t)− Φs,n−1(t)2∆x(3.34)∂Φe(x, t)∂x|xn ≈Φe,n+1(t)− Φe,n−1(t)2∆x(3.35)∂ ln ce(x, t)∂x|xn ≈ln ce,n+1(t)− ln ce,n−1(t)2∆x(3.36)using a central differencing scheme. Together, the temperature equation inthe electrodes can be discretized as,ρiCp,i∂T (x, t)∂t= λiTn−1 − 2Tn + Tn + 1∆x+Qohm,n +Qrxn,n +Qrev,n(3.37)3.2.2 Implementation of Boundary and InterfaceConditionsBoundary conditions require certain variables being evaluated at theends of segments. For example, consider the electrolyte potential Φe at the423.2. Numerical Implementationinterface between the anode and the negative current collector,Φe(x, t)∣∣∣x=xn= 0 (3.38)In order to recover such value, the ghost point technique is used, as shownin Fig. 3.4.Figure 3.4: Interpolation technique to recover interface valuesThe discretized equation is thus,Φe,N + Φe,N+12= 0 (3.39)We can apply the same approach to continuity and interface conditions.Consider the electrolyte concentration ce. Since electrolyte is present in allof cathode, separator, and anode, continuity of the solution for both theconcentration ce and the potential Φe have to be enforced at the cathode-separator and the separator-anode junctions. The easiest way would beto use the ghost point technique. For example, at the cathode-separator433.2. Numerical Implementationjunction, the continuity condition for the electrolyte concentration is,ce(x, t)∣∣∣x=xˆ−p= ce(x, t)∣∣∣x=xˆ+p(3.40)which can be discretized as,ce,pN + ce,pN+12=ce,s0 + ce,s12(3.41)where suffixes pN and pN+1 represent the last and the ghost point of the cath-ode, and s1 and s0 represent the first and the ghost point of the separator.Fig. 3.5 is a pictorial description of this interface.Figure 3.5: Electrolyte continuity across the cathode and the separatorSimilarly, continuity of fluxes across interfaces is also enforced. Consider theequation of diffusion coefficient Deff at the same cathode-separator junction,−Deff,p∂ce(x, t)∂x∣∣∣x=xˆ−p= −Deff,s∂ce(x, t)∂x∣∣∣x=xˆ+p(3.42)which can be discretized as,− Deff,pN +Deff,pN+12ce,pN+1 − ce,pN∆x= −Deff,s0 +Deff,s12ce,s1 − ce,s0∆x(3.43)Again, suffixes pN and pN+1 represent the last and the ghost point of the443.2. Numerical Implementationcathode, and s1 and s0 represent the first and the ghost point of the sepa-rator. Notice that all boundary conditions discretized using the ghost pointtechnique retains second-order accuracy.3.2.3 Time Discretization and Newton’s MethodWith spatial discretization completed, we now proceed to discretize time.Backward Euler (BE) time-stepping is chosen in order to maintain stabilityat each time-step while keeping an approriate size of ∆t. The example belowshows a fully discretized solid particle concentration equation,ck+1s,n − cks,n∆t=Dspr2mr2m+1/2(ck+1s,n,m+1 − ck+1s,n,m)− r2m−1/2(ck+1s,n,m − ck+1s,n,m−1)∆r2(3.44)where k indicates the current time-step and k = 0 refers to the initial con-dition.With every of the P2D equations discretized in space and time, we canarrange all the variables into a vector uk+1,uk+1 =[ck+1e,p ck+1e,s ck+1e,n ck+1s,p ck+1s,n jk+1p jk+1n Φk+1s,p Φk+1s,nΦk+1e,p Φk+1e,s Φk+1e,n Tk+1a Tk+1p Tk+1s Tk+1n Tk+1z]Twhere k represents time-step. The goal is that, with known initial conditionu0, we want to find uk+1 for k = 0, 1, 2, ...,K, which is described by thematrix equation,A · uk+1 + v − ucur = 0 (3.45)453.2. Numerical Implementationwhere A is a constant coefficient matrix that takes into account all linearpart of every P2D equation, v is a vector consisting all nonlinear part, anducur is a vector containing information about the current time-step,ucur =[cke,p cke,s cke,n cks,p cks,n 0 0 0 0 0 0 0T ka Tkp Tks Tkn Tkz]TSince the majority of the P2D equations are nonlinear, matrix Newton’smethod is used to find the root. To use Newton’s method, letF = A · uk+1 + v − ucur (3.46)We compute matrix J , the derivative of F with respect to uk+1,J = A+Dv (3.47)where Dv is the Jacobian matrix of vector v, namely, Dvij =∂vi∂uj. Approx-imation of the P2D equations is now possible with the following algorithm,Algorithm: Backward Euler Approximation to P2D equationsconstruct matrix Aconstruct vector ucurfor time-steps k = 0, 1, 2, ..., Tf doconstruct vector v and matrix Dvconstruct matrix Jcalculate residual=A · uk+1 + v − ucurwhile residual < tol do463.2. Numerical Implementationfind solution y to J · y = residualupdate uk+1 = uk+1 − yupdate vector v and matrix Dvupdate matrix Jupdate residual=A · uk+1 + v − ucurend whileupdate ucur with values of uk+1end forThe residual here is calculated using the formula,residual = ||A · uk+1 + v − ucur||∞ (3.48)where we are using the maximum norm. Numerically, Newton iteration willcontinue until the residual is smaller than a specified tolerance tol. Notethat Newton iteration can fail to converge or find a different root to the onesought after if the function has many inflection points or if the initial guessis not close enough. For our purposes, we will not run into these situationsas our initial guess computed explicitly using the initial condition given inTable II will be sufficiently close in our BE solution.3.2.4 Implementation ResultsSimulation results were obtained using MATLAB R2018a on a Windows10@1.8GHz PC with 16GB of RAM for the experimental battery parametersin Table IV with a cutoff voltage of 2.5V and environmental temperature473.2. Numerical Implementationof 298.15K. For the proposed chemistry, the 1C value is ≈ 30 A/m2. Thedefault discretization sets ∆x = 1×10−6 and ∆r = 0.5×10−6 unless specifiedotherwise. The battery voltage is calculated by taking the difference betweenthe solid partial potential of the first segment of the cathode and that of thelast segment of the anode,V = Φs,p1 − Φs,nNIn the first scenario shown in Fig. 3.6, a 1C discharge simulation witha fixed value of h = 1W/(m2K) is performed and the average temperaturefor each section is plotted. Since the thermal conductivity coefficients areextremely high (λa = 237,λp = 2.1, λs = 0.16, λn = 1.7, and λz = 401) withrespect to the length scale of each section (O(10−5)), heat diffusion is suffi-ciently fast through the entire battery so that the temperature is virtuallythe same across all five sections. Therefore, in subsequent discussions, thebattery temperature simply refers to the average temperature across all fivesections.In the second scenario shown in Fig. 3.7, 1C discharge simulations arecompared for a wide range of heat exchange coefficient h. As expected,decreasing the value of the heat exchange coefficient h leads to a more insu-lated battery and thus a faster increase of the cell temperature. Moreover,due to the coupling of all of the governing equations, it is possible to notethe influence of different temperatures on the cell voltage.In the third scenario shown in Fig. 3.8, for a fixed value of h = 1W/(m2K),different discharge cycles are compared at 1C, 2C, and 5C. According to the483.2. Numerical ImplementationFigure 3.6: Average temperature of each battery section in 1C dischargewith h = 1W/(m2K)different applied currents, the temperature rises in different ways. It is in-teresting to note the high slope of the temperature during the 5C discharge,mainly due to the electrolyte concentration ce being driven to zero in thepositive electrolyte by the high discharge rate.493.2. Numerical ImplementationFigure 3.7: 1C discharge cycle run under different heat exchange coefficients:h = 0.01, h = 1, and h = 1003.2.5 Verification of Implementation ResultsFrom Taylor series, we observe that our discretization of time using BEis first-order accurate,∂f∂t∣∣∣tk+1≈ fk+1 − fk∆t+O(∆t) (3.49)503.2. Numerical ImplementationFigure 3.8: Discharge cycle run under 1C, 2C, and 5Cand our discretization of space for both first and second derivatives usingcentral-difference is second-order accurate,∂f∂x≈ fn+1 − fn2∆x+O(∆x2) (3.50)∂2f∂x2≈ fn−1 − fn + fn+1∆x2+O(∆x2) (3.51)We can now verify the implementation results by explicitly calculating theorder of convergence.513.2. Numerical ImplementationSuppose we had an exact solution uexact to the P2D equations and letuapprox denote our approximation to the exact solution. Without the loss ofgenerality, the global error E is calculated by the absolute value of the dif-ference between the cell voltage vs. time profile from the exact solution andthat from the approximate solution, averaged over the length of simulationtime,E =∫ tfti||Vexact(t)− Vapprox(t)||dttf − ti (3.52)Since we have O(∆t) error from Backward Euler and O(∆x2) error fromfinite difference, the order of the total error in our implementation is,E = O(∆t) +O(∆x2) (3.53)Our method is said to converge if this error vanishes as ∆x→ 0 and ∆t→ 0.We now use our simulation results to illustrate that our implementationis indeed first-order in time and second-order in space. To confirm that it isfirst-order in time, we fix the number of grid points in both x and r, so thatthe spatial error O(∆x2) remains constant and will not affect convergence.We then divide the time-step size ∆t by 2 successively. The error from time-stepping Etime is calculated by taking the maximum norm of each successiveapproximation normalized over total simulation time,Etime =∫ tfti||V∆t/2 − V∆tdt||tf − ti (3.54)523.2. Numerical Implementationand the estimated convergence rate can then be written as,CR = log2(∫ tfti||V∆t/4 − V∆t/2dt||∫ tfti||V∆t/2 − V∆tdt||)(3.55)with CR = 1 for first-order and CR = 2 for second-order. First-order con-vergence in time is clearly seen in Table 3.1∆t Etime CR201.1959× 10−21.0325105.8463× 10−31.019252.8845× 10−32.5Table 3.1: First-order convergence in time for Backward Euler methodThe situation is more complicated when comes to checking the conver-gence rate in space, as there are two dimensions x and r with ∆x 6= ∆r.Therefore, we will need to check each one individually. First, we fix ∆tand ∆r so that convergence is dependent solely on ∆x. We then half ∆xsuccessively. The error Ex from discretizing x is calculated by taking themaximum norm of each successive approximation normalized over total sim-ulation time,Ex =∫ tfti||V∆x/2 − V∆xdt||tf − ti (3.56)and the estimated convergence rate can then be written as,CR = log2(∫ tfti||V∆x/4 − V∆x/2dt||∫ tfti||V∆x/2 − V∆xdt||)(3.57)533.2. Numerical Implementation∆x Ex CR1× 10−66.1631× 10−71.36190.5× 10−62.3979× 10−70.25× 10−6Table 3.2: Error and convergence for discretization in xThe results are shown in Table 3.2. Different from our expectation, the CRhere is less than 2, which may be due to the extremely small magnitude oferror that interferes with MATLAB’s intrinsic machine precision.Lastly, we fix ∆t and ∆x so that convergence is dependent solely on ∆r.We then half ∆r successively. Simiarly, the error Er from discretizing r iscalculated by taking the maximum norm of each successive approximationover total simulation time,Er =∫ tfti||V∆r/2 − V∆rdt||tf − ti (3.58)and the estimated convergence rate can then be written as,CR = log2(∫ tfti||V∆r/4 − V∆r/2dt||∫ tfti||V∆r/2 − V∆rdt||)(3.59)The results are shown in Table 3.3. Here, CR is asymptotically approachingthe expected second-order convergence rate as we half ∆r. Notice that, of∆t, ∆x, and ∆r, it is ∆r which has the most impact on the accuracy of ourimplementation.543.3. Model Reduction∆r Er CR0.5× 10−61.8701× 10−23.05940.25× 10−62.2434× 10−32.09600.125× 10−65.2474× 10−40.0625× 10−6Table 3.3: Error and convergence for discretization in r3.3 Model ReductionThe ultimate goal of P2D model simulation is to implement it on ad-vanced battery management systems (ABMS). ABMS anticipate problemsthrough online fault diagnosis which can prevent damage, ensure safety, min-imize charging time, and slow down battery aging. These are possible onlyif model simulations are extremely fast. The full P2D model is currentlythe most accurate electrochemical model but the computational time is toolong to be implemented on ABMS. To achieve a better trade-off betweenaccuracy and computational time, four different model reductions are pro-posed: the two-parameters approximation model, the temperature-reductionmodel, the Φs-reduction model, and the mixed-reduction model. The idealreduced model can achieve less than 1% error when compared with the fullmodel at any discharge rate. The convergence of each approximate modelis assessed by comparing its cell potential vs. time profiles under differentapplied current rates Iapp(t) with that of the full P2D model, as shown inFig. 3.9.553.3. Model ReductionFigure 3.9: Comparison of the three reduced models with the full model. a)1C rate comparison. b) 2C rate comparison. c) 5C rate comparison.The improved efficiency of the four reduced models is demonstrated inTable 3.4 by comparing their individual effective simulation time with thatof the full model.563.3. Model ReductionSimulation Duration Effective Simulation TimeFull model 3590s 341.016sTwo-parameters approximation model 3540s 203.716sReduced temperature model 3580s 161.194sReduced Φs model 3600s 261.079sMixed reduction model 3600s 123.921sTable 3.4: Timing comparisons of different models running under 1C, with∆t = 10, ∆x = 1× 10−6, ∆r = 0.5× 10−6, and h = 13.3.1 Two-Parameters Approximation ModelModel overviewRecall that diffusion inside solid spherical particles is described by Fick’slaw,∂cs(r, t)∂t=1r2∂∂r[r2Dsp∂cs(r, t)∂r](3.60)with boundary conditions,∂cs(r, t)∂r∣∣∣r=0= 0,∂cs(r, t)∂r∣∣∣r=Rp= −j(x, t)Dseff(3.61)In this model, a major source of computational burden comes from thepseudo-second-dimension (r).In the two-parameters approximation model, concentration profiles in-side the particle are assumed to be quadratic in r,cs(r, t) = a(t) + b(t)r2R2p(3.62)where a(t) and b(t) are to be determined. Substituting Eq. 3.62 into Eq.573.3. Model Reduction3.60, we obtain,∂a(t)∂t+r2R2p∂b(t)∂t− 6Dseffb(t)R2p(3.63)The boundary condition at r = 0 is automatically satisfied. The boundarycondition at r = Rp becomes,2DseffRpb(t) = −j (3.64)We are interested in the average concentration and surface concentration.Thus, a(t) and b(t) are expressed in terms of volume-average concentrationcavgs (x, t) and surface concentration c∗s(x, t). The volume-averaged concen-tration is given by,cavgs (x, t) =∫ Rpr=03r2R2pcs(r, t)d( rRp)(3.65)Applying Eq. 3.62 in Eq. 3.65, we get,cavgs (x, t) = a(t) +35b(t) (3.66)The surface concentration c∗s(x, t) is obtained by substituting r = Rp in Eq.3.62,c∗s(x, t) = a(t) + b(t) (3.67)583.3. Model ReductionEq. 3.66 and 3.67 are solved to obtain,a(t) = −32c∗s(x, t) +52cavgs (x, t) (3.68)b(t) = −52cavgs (x, t) +52c∗s(x, t) (3.69)Now we can substitute a(t) and b(t) into Eq. 3.62,cs(r, t) = −32c∗s(x, t) +52cavgs (x, t) +(− 52cavgs (x, t) +52c∗s(x, t)) r2R2p(3.70)We need two equations to evaluate the average concentration cavgs (x, t) andthe surface concentration c∗s(x, t). The former can be evaluated by volumeaveraging the entire governing Eq. 3.60,∫ Rpr=03r2R2p[∂cs∂t−Dseff1r2∂∂r(r2∂cs∂r)]d( rRp)= 0 (3.71)Substituting Eq. 7 and 13, and evaluating, we have,∂cavgs (x, t)∂t= −3j(x, t)Rp(3.72)The surface concentration is obtained by evaluating the boundary conditionat r = Rp. Applying Eq. 3.92 to Eq. 3.62, we have,c∗s(x, t)− cavgs (x, t) = −RpDspj(x, t)5(3.73)In summary, the original diffusion inside solid-particles is now approxi-593.3. Model Reductionmated by means of average and surface concentration,∂cavgs (x, t)∂t=− 3j(x, t)Rpc∗s(x, t)− cavgs (x, t) =−RpDspj(x, t)5This reduction leads to a one-dimensional problem in x by removing thepseudo-second-dimension r. In terms of computation, 47.63% of variablesare eliminated for any specified ∆x and ∆r.Results analysisFrom Fig. 3.9 we observe that for medium (2C) and high (5C) dischargerates, the model simulation ends prematurely, primarily because the elec-trolyte concentration ce being driven to zero in the positive electrode by thehigh discharge rate.The performance of the model is evaluated by the absolute difference ofcell potential profile averaged over total simulation time with the full model,Etp =∫ tfti||Vreduced − Vfulldt||tf − ti (3.74)Note that since the reduced model ends prematurely, tf is normalized byextending the earlier ending model with V = 0 until it has the same lengthof time as the full model. The performance of the two-parameters approxi-mation model can be seen in Table 3.5,In terms of accuracy, the two-parameters approximation model meetsour criterion of less than 1% error only for discharge rate < 1C. In terms603.3. Model ReductionEtp1C 0.00884582C 0.0583245C 0.28877Table 3.5: Performance of the two-parameters approximation modelof effective simulation time (Table 3.4), this model improves computationalefficiency by 40%.3.3.2 Temperature-Reduction ModelModel overviewFrom Fig. 3.6, the temperature is shown to be constant in all sections,in other words T is constant in x. Thus, we can reduce computation time bysimply having one single global temperature variable T instead of a differentT variable on each grid point. This model reduction eliminates 8.99% ofvariables for any specified ∆x and ∆r.Results analysisSimilar to the two-parameters approximation, the temperature-reductionmodel ends prematurely at 2C and 5C, as seen in Fig. 3.9. However at 5C,its potential vs. time profile is clearly closer to the full model. To confirmmodel performance, we use the same metric,Etr =∫ tfti||Vreduced − Vfulldt||tf − ti (3.75)613.3. Model ReductionAgain, tf is normalized by extending the earlier ending model with V = 0until it has the same length of time as the full model. The performance ofthe temperature-reduction model can be seen in Table 3.6,Etr1C 0.00220762C 0.0515455C 0.16308Table 3.6: Performance of the temperature-reduction modelIn terms of accuracy, the temperature-reduction model meets our cri-terion of less than 1% error only for discharge rate < 1C. In terms ofeffective simulation time (Table 3.4), although this model eliminates lessvariables than the previous reduced model, it achieves the highest compu-tational efficiency. This is because temperature is embeded in every of theP2D equations. Reducing all temperature into one single variable convertsmany vector and matrix operations into scalar calculations, and thus greatlyshortens simulation time.3.3.3 Φs-Reduction ModelModel overviewThe effective diffusivity of electrolyte and solid particles, Deff and Dseff,are of the same magnitude. However, because the length of electrolyte(8× 10−5m for cathode and 8.8× 10−5m for anode) is ∼ 20 times of that ofthe solid particles (4× 10−6m in diameter), actual diffusion is much fasterin the latter. Consequently, the solid concentration cs is nearly uniform in623.3. Model Reductionboth cathode and anode, and thus there is little solid potential (Φs) varia-tion in each section. Therefore, instead of having Np and Nn identical valuesof Φs in the cathode and anode, respectively, we can reduce the model tohave only two Φs, one for each section. This eliminates 6.75% of variablesfor any specified ∆x and ∆r.Results analysisFrom Fig. 3.9, we observe that the Φs-reduction model overlaps the fullmodel almost completely for all three discharge rates. To confirm modelperformance, we again use the standard,Espr =∫ tfti||Vreduced − Vfulldt||tf − ti (3.76)where tf is normalized similar to the previous two reduced models by ex-tending the earlier ending model with V = 0 until it has the same length oftime as the full model. The performance of the Φs-reduction model can beseen in Table 3.7,Espr1C 0.00140622C 0.00241865C 0.004695Table 3.7: Performance of the Φs-reduction modelIn terms of accuracy, the Φs-reduction model exceeds our criterion of lessthan 1% error for all of low (1C), medium (2C), and high (5C) dischargerates. In terms of effective simulation time (Table 3.4), this model improves633.3. Model Reductioncomputational efficiency by 23%.3.3.4 Mixed-Reduction ModelModel overviewLastly we try to combine all model reductions into a single model. Thatis, the mixed-reduction model incorporates the above three model reduc-tions: the two parameter approximation model, the temperature-reductionmodel, and the Φs-reduction model. In terms of computation, 63% of vari-ables are eliminated for any specified ∆x and ∆r.Results analysisFrom Fig. 3.9 we observe that the mixed-reduction model simulationends much prematurely and produces results with the largest offset whencompared to the full model. The performance deteriorates quickly at highercurrent rates.To confirm model performance, we again use the standard,Emr =∫ tfti||Vreduced − Vfulldt||tf − ti (3.77)Same as before, tf is normalized by extending the earlier ending modelwith V = 0 until it has the same length of time as the full model. Theperformance of the mixed-reduction model can be seen in Table 3.,In terms of accuracy, the mixed-reduction model meets our criterion ofless than 1% error only for discharge rate < 1C. At 2C and 5C, this modelyields errors of 10% and 39%. This large error is expected as this model643.4. ConclusionEmr1C 0.00926872C 0.101255C 0.3893Table 3.8: Performance of the mixed-reduction modeleliminates a high percentage of fundamental variables. In terms of effectivesimulation time (Table 3.4), this model improves computational effiency by64%, the largest efficiency increase of all. Applications that require only lowcurrent rates but fast response time on battery monitoring may incorporatethis model into their BMS.3.4 ConclusionThis chapter describes a detailed procedure for the numerical imple-mentation of the pseudo-two-dimensional (P2D) model. The treatment ofboundary conditions is addressed with particular attention to the inter-face conditions across the different sections of the battery. The simulationsdemonstrate high numerical stability for different operating scenarios.Four model reductions with elimination of different variables to reducecomputational complexity are implemented and tested. All reduced modelshave effective simulation time shorter than that of the full model, withmixed-reduction model being the fastest, followed by reduced-temperaturemodel, two-parameters approximation model, and lastly reduced Φs model.The results of this work show the promise of the proposed framework asa reliable and efficient MATLAB-based program for the P2D model simula-653.4. Conclusiontion. Further developments such as code migration into a compiled language(e.g. C#, Java) can only improve the current performance. As the proposedsimulations are written in standard serial mode, the computation time couldbe reduced by at least a factor of ten by using a multicore CPU using par-allel DAE solvers. Modern versions of MATLAB have easy-to-implementatbuilt-in options for distributing calculations among multiple cores on a singleCPU, and among multiple CPUs.66Chapter 4Summary and Future WorkHere we provide a summary for each of the two approaches taken anddiscuss future extensions to the work done.4.1 First Approach: Physical ExperimentsIn this approach, lithium-ion batteries are cycled under different cur-rent rates, while their temperature, voltage, current, and cell capacity aremeasured. Battery aging, which is accelerated by high current rates, is illus-trated with plotting and analysis. Two predictive models that can estimatelithium-ion battery state-of-health are built from the data generated. Fu-ture work that may be interesting includes distinguishing battery aging dueto cycling (discharge/charge) from purely temperature variation. This canbe done by setting up two treatments: the first treatment involves cyclingbatteries with different current rates while recording their temperature, andthe second treatment involves no discharge/charge but thermocycles bat-teries using the temperature profile generated from the first treatment. Anecessary equipment needed for this work is heating tape with a customprogrammed controller.674.2. Second Approach: Mathematical Modelling4.2 Second Approach: Mathematical ModellingIn this approach, we provided a robust implementation of the pseudo-two-dimensional (P2D) model and discussed four potential model reduc-tions that greatly shorten simulation time. Each model reduction achievesa different standard of accuracy, which can be programmed onto batterymangement system of different applications. An interesting future workthat may further improve computational efficiency includes developing asplit-step solver. In this split-step solver, with a and c given by the initialcondition, first update the nonlinear parts j, s, Φ, and φ with Newton’smethod, then update T and a explicitly using the forward Euler method.Lastly, with j and T known from the previous two steps, update c using theIMEX (implicit-explicit) method. Overall, this split-step solver tailors themost efficient method to each equation: implicit method for the nonlinearequations and explcit method for the linear equations. The simulation timeis thus expected to further shorten.68Bibliography[1] Atebamba, J.M. et. al., On the Interpretation of Measured ImpedanceSpectra of Insertion Cathodes for Lithium-Ion Batteries Journal of TheElectrochemical Society. 157, 11 (2010), 1218-1228.[2] Cai, L. and White, R.E., Mathematical Modeling of a Lithium-ion Bat-tery with Thermal Effects in COMSOL Inc. Multiphysics (MP) SoftwareJournal of Power Sources. 196, 14 (2011), 5985 - 5989.[3] Dong, H. et. al., Lithium-ion Battery SOH Monitoring and RemainingUseful Life Prediction based on Support Vector Regression-Particle Fil-ter Journal of Power Sources. 271 (2014), 114-123.[4] Dong H.J., Numerical Modeling of Lithium Ion Battery for Predict-ing Thermal Behavior in a Cylindrical Cell Current Applied Physics.14 (2014), 196-205.[5] Doyle C.M., Design and Simulation of Lithium Rechargeable BatteriesPh.D. Dissertation, University of California, Berkeley. August 1995[6] Fernandez, I.J. et. al., Capacity Fade and Aging Models for ElectricBatteries and Optimal Charging Strategy for Electric Vehicles. Energy.60 (2013), 35-43.69Bibliography[7] Hausbrand, R. et. al., Fundamental Degradation Mechanisms of LayeredOxide Li-ion Battery Cathode Materials - Methodology, Insights andNovel Approaches. Materials Science & Engineering B. 192 (2015), 3-25.[8] Jokar A. et. al., Review of Simplified Pseudo-two-Dimensional Modelsof Lithium-ion Batteries Journal of Power Sources. 327 (2016), 44-55.[9] Kemper P. et. al., Simplification of Pseudo two Dimensional BatteryModel using Dynamic Profile of Lithium Concentration Journal of PowerSources. 286 (2015), 510-525.[10] Kespe M. and Nirschl H., Numerical Simulation of Lithium-ion BatteryPerformance Considering Electrode Microstructure International Jour-nal of Energy Research. 39 (2015), 2062-2074.[11] Love, C. et. al., State-of-Health Monitoring of 18650 4S Packs witha Single-Point Impedance Diagnostic. Journal of Power Sources. 266(2014), 512-519.[12] Nagaoka N., A Numerical Model of Lithium-ion Battery for a LifeEstimation International Universities’ Power Engineering Conference(UPEC). (2013), 1-6.[13] Omar, N. et. al., Lithium Iron Phosphate Based Battery - Assessmentof the Aging Parameters and Development of Cycle Life Model AppliedEnergy. 113 (2014), 1575-1585.[14] Ramadesigan, V. et. al., Efficient Reformulation of Solid-Phase Diffu-70sion in Physics-based Lithium-ion Battery Models Journal of The Elec-trochemical Society. 157, 7 (2010), A854.[15] Salkind A.J. et. al., Determination of State-of-Charge and State-of-Health of Batteries by Fuzzy Logic Methodology Journal of PowerSources. 80 (1999), 293-300.[16] Samad N.A. et. al., Battery Capacity Fading Estimation Using a Force-Based Incremental Capacity Analysis Journal of Electrochemical Society.163, 8 (2016), 1584-1594.[17] Stroe D.I. et. al., Diagnosis of Lithium-ion Batteries State-of-HealthBased on Electrochemical Impedance Spectroscopy Technique IEEE En-ergy Conversion Congress and Exposition (ECCE). (2014), 4576-4582.[18] Torchio M. et. al., LIONSIMBA: A Matlab Framework Based on aFinite Volume Model Suitable for Li-Ion Battery Design, Simulation,and Control Journal of Electrochemical Society. 163 (2016), 1192-1205.[19] Williard N. et. al., Comparative Analysis of Features for DeterminingState-of-Health in Lithium-ion Batteries International Journal of Prog-nostics and Health Management. 4, 1 (2013), 14-20.[20] Zhang W. et. al., A Review of the Electrochemical Performance of Al-loy Anodes for Lithium-ion Batteries Journal of Power Sources. 196, 1(2013), 13-24.7172Appendix A. AppendixAppendix AAppendixTable I. Additional equationsOpen circuit potential (thermal dependence)Up = Up,ref + (T (x, t)− Tref)∂Up∂T |TrefUn = Un,ref + (T (x, t)− Tref)∂Un∂T |TrefEntropy change∂Up∂T∣∣∣Tref= −0.001(0.199521039−0.928373822θp+1.364550689000003θ2p1−5.661479886999997θp+11.47636191θ3p−0.6115448939999998θ3p+3.048755063θ4p)∂Un∂T∣∣∣Tref=(0.001(0.005269056+3.299265709θn−91.79325798θ2n+1004.911009θ3n(1−48.09287227θn+1017.234804θ2n−5812.278127θ4n+19329.7549θ5n−37147.8947θ6n+38379.18127θ7n−10481.80419θ3n+59431.3θ4n−195881.6488θ5n+374577.3152θ6n−16515.05308θ8n)−385821.1607θ7n+165705.8597θ8n))Open circuit potential (reference value)Up,ref =−4.656+88.669θ2p−401.119θ4p+342.909θ6p−462.471θ8p+433.434θ10p−1+18.933θ2p−79.532θ4p+37.311θ6p−73.083θ8p+95.96θ10pUn,ref =(0.7222 + 0.1387θn + 0.029θ0.5n − 0.0172θn + 0.0019θ1.5n+0.2808e0.9−15θn − 0.7984e0.4465θn−0.4108)θp =c∗s,p(x,t)cmaxs,pθn =c∗s,n(x,t)cmaxs,nHeat source terms (anode and cathode)Qohm = σeff,i(∂Φs(x,t)∂x)2+ κeff,i(∂Φe(x,t)∂x)2+2κeff,iRT (x,t)F (1− t+)∂ ln ce(x,t)∂x ∂Φe(x,t)∂xQrxn = Faij(x, t)ηi(x, t)Qrev = Faij(x, t)T (x, t)∂Ui∂T∣∣∣Tref73Appendix A. AppendixTable I. Additional equations (continued)Heat source terms (separator)Qohm = κeff,i(∂Φe(x,t)∂x)2+2κeff,iRT (x,t)F (1− t+)∂ ln ce(x,t)∂x ∂Φe(x,t)∂xVarious coefficientsDeff,i = bruggii × 10−4 × 10−4.43− 54T (x,t)−229−5×10−3ce(x,t)−0.22×10−3ce(x,t)κeff,i = bruggii × 10−4 × ce(x, t)(− 10.5 + 0.668× 10−3ce(x, t)+0.494× 10−6c2e(x, t) + (0.074− 1.78× 10−5ce(x, t)−8.86× 10−10c2e(x, t))T (x, t) + (−6.96× 10−5+2.8× 10−8ce(x, t))T 2(x, t))2keff,i = kie−EkiaR(1T (x,t)− 1Tref)Dseff,i = Dsi e−EDsiaR(1T (x,t)− 1Tref)σeff,i = σi(1− i − f,i)Υ := 2(1−t+)RF74Appendix A. AppendixTable II. List of parameters used in simulationParameter Descriptioncinite [mol/m3] Initial concentration in the electrolytecavg,inits [mol/m3] Initial solid-phase concentrationcmaxs [mol/m3] Maximum solid-phase concentrationDi [m2/s] Electrolyte diffusivityDsi [m2/s] Solid-phase diffusivityki [m2.5/(mol0.5s)] Reaction rate constantli [m] ThicknessRp,i [m] Particle radiusρi [kg/m3] DensityCp,i [J/(kg·K)] Specific heatλi [W/(m·K)] Thermal conductivityσi [S/m] Solid-phase conductivityi [-] Porosityai [m2/m3] Particle surface area to volumeEDsia [J/mol] Solid-phase diffusion activation energyEkia [J/mol] Reaction constant activation energybrugg [-] Bruggeman’s coefficientF [C/mol] Faraday’s constantR [J/(mol·K)] Universal gas constantt+ [-] Transference numberf,i [-] Filler fraction75Appendix A. AppendixTable II. List of parameters used in simulation (continued)Parameter Aluminum CC Cathode Separatorcinite - 1000 1000cavg,inits - 25751 -cmaxs - 51554 -Di - 7.5× 10−10 7.5× 10−10Dsi - 10−14 -ki - 2.334× 10−11 -li 10−5 8× 10−5 2.5× 10−5Rp,i - 2× 10−6 -ρi 2700 2500 1100Cp,i 897 700 700λi 237 2.1 0.16σi 3.55× 107 100 -i - 0.385 0.724ai - 885,000 -EDsia - 5000 -Ekia - 5000 -brugg - 4 4F 96485 96485 96485R 8.314472 8.314472 8.314472t+ 0.364 0.364 0.364f,i - 0.025 -76Appendix A. AppendixTable II. List of parameters used in simulation (continued)Parameter Anode Carbon CCcinite 1000 -cavg,inits 26128 -cmaxs 30555 -Di 7.5× 10−10 -Dsi 3.9× 10−14 -ki 5.031× 10−11 -li 8.8× 10−5 10−5Rp,i 2× 10−6 -ρi 2500 8940Cp,i 700 385λi 1.7 401σi 100 5.96× 107i 0.485 -ai 723,600 -EDsia 5000 -Ekia 5000 -brugg 4 -F 96485 96485R 8.314472 8.314472t+ 0.364 0.364f,i 0.0326 -77


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items