energies Article Multiphysics Modeling for Detailed Analysis of Multi-Layer Lithium-Ion Pouch Cells NanLin1,2 ,FridolinRöder1,2andUlrikeKrewer1,2,* 1 MechanicalEngineeringDepartment,InstituteofEnergyandProcessSystemsEngineering, TechnischeUniversitätBraunschweig,Franz-Liszt-Str.35,D-38106Braunschweig,Germany; [email protected](N.L.);[email protected](F.R.) 2 BatteryLaboratoryBraunschweig,TUBraunschweig,D-38106Braunschweig,Germany * Correspondence:[email protected];Tel.:+49-531-391-3030 (cid:1)(cid:2)(cid:3)(cid:1)(cid:4)(cid:5)(cid:6)(cid:7)(cid:8)(cid:1) (cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7) Received:14October2018;Accepted:29October2018;Published:1November2018 Abstract: Multiphysicsmodelingpermitsadetailedinvestigationofcomplexphysicalinteractions andheterogeneousperformanceinmultipleelectro-activelayersofalarge-formatLi-ioncell. Forthis purpose, a novel 3D multiphysics model with high computational efficiency was developed to investigate detailed multiphysics heterogeneity in different layers of a large-format pouch cell at various discharge rates. This model has spatial distribution and temporal evolution of local electric current density, solid lithium concentration and temperature distributions in different electro-activelayers,basedonarealpouchcellgeometry. Otherthanpreviousmodels,weresolvethe dischargeprocessesatvariousdischargeC-rates,analyzinginternalinhomogeneitybasedonmultiple electro-activelayersofalarge-formatpouchcell. Theresultsrevealthatthestronginhomogeneityin multiplelayersatahighC-rateiscausedbythelargeheatgenerationandpoorheatdissipationin thedirectionthroughthecellthickness. Thethermalinhomogeneityalsostronglyinteractswiththe localelectrochemicalandelectricperformanceintheinvestigatedcell. Keywords: 3Dmultiphysicsmodel;lithium-ionbattery;batterydesign;heterogeneity 1. Introduction Thereisagrowingdemandforadvancedenergystoragesystems.Lithium-ionbatteries(LiBs)have becomeahigh-demandenergytechnologyduetotheirhighenergydensity.Asaconsequence,theyare widelyappliedinenergystorageandsupplyrangingfromelectricalvehicles(EVs),hybridelectrical vehicles(HEVs),andsmall-sizeddevicesinsatellites[1–3]. However,theneedtoovercomeexisting barriers,suchasreducingcostandfurtherincreasingenergydensity,motivatestheinvestigationof the impact of macroscopic design features on internal physical processes and their interaction [4]. Thoseaspectsareparticularlyimportantinlarger-formatcells,e.g.,pouchcellsorcylindricalcells. Such large-format cells are made of multiple layers to achieve high energy density. Therefore, forimprovementoflarge-formatcells,itisnecessarytodeeplyinvestigateasinglemulti-layercell, asadetailedinsightofheterogeneousmultiphysicsprocessesinthecellishelpfulforadirecteddesign andoptimizationoflarge-formatcells. Mathematicalmodelscanbeusedtodescribevariouscoupledphysicalprocessestoevaluatetheir interaction. SuchmultiphysicsmodelsofLiBscanbeappliedfordetailedanalysisofinternalphysical processesandbatteryperformance. Higherdimensionalmultiphysicalmodelscanbeusedtoreveal thoseaspectsinlarge-formatcellsandevaluatelocalheterogeneityfordifferentoperatingconditions, e.g.,chargeandenvironmentaltemperature. To establish multiphysics modeling of LiBs, prior studies applied different approaches and coupling methodologies. Several fully coupled 2D electrochemical-thermal models have been Energies2018,11,2998;doi:10.3390/en11112998 www.mdpi.com/journal/energies Energies2018,11,2998 2of26 developed to investigate multiphysics processes in a specialcross-sectionedplane. They are used to analyze the electrochemical effects of current collectors and tabs [5], temperature and current distributionsacrossthebatterysurface[6,7]andthepseudo-3Dporouselectrodemodelfordetailed ionictransport[8]. Duetothelimitationofmodeldimensions,thesemodelscannotentirelysimulate the real internal physical processes with local heterogeneity. 3D multiphysics modeling of LiBs, with different coupling methods, has been employed to investigate cell multiphysics with 3D cell structures[9–11].However,thesepapersdidnotdiscussabouttheinfluencesofelectro-activelayer numbersandstructuresontheheterogeneousperformanceoflarge-formatpouchcells. Thediscrete structuresinlarge-formatcellshaveastrongeffectonbatteryperformance. Fortheinvestigationofbatterypacksormodules,Guoetal. [12]developeda3Dmultiphysics model. Thismodelhasbeenusedtostudythethermalandelectrochemicalperformanceinamodule madeofthreepouchcells. Despitedetailedillustrationofcontourresultsofelectricpotentialand temperature,themodeldidnotshowtheinternalinhomogeneityofeachsinglecell. Pannalaetal. [13] alsopresenteda3Dbatterypackmodel, whichwasusedtosimulateallbatterycomponentsfora detailedevaluationofthethermalandmechanicalsafetyriskunderadverseconditions. Although multiphysics modeling of large-format cells has already enabled a detailed investigation of cell components and the impacts of the battery materials on its performance, detailed analysis of inhomogeneousperformancesinasinglelarge-formatpouchcellwithmultipleelectro-activelayers wasbarelyaddressed. Moreover,toreducecomputationalcost,internalmaterialhomogeneityisoften assumed. Therefore,recentmultiphysicsmodelsprovidefewresultsforinternalheterogeneitybased onelectro-activelayerstructuresofasinglelarge-formatcell,includingdynamicdescriptionsofsuch heterogeneousmultiphysicsprocesses. Bycombiningtwohierarchicalframeworks[14,15]forgoodmanagementofcouplinginterfaces andsubmodelsolvers,a3Dmultiphysicsmodelisdevelopedforthepredictionandinvestigationof internalinhomogeneityregardingcellstructures,whichisparticularlyapplicabletolarge-formatpouch cellswithmultipleelectro-activelayers. Ithasbeenappliedforthesimulationofa12Ahpouchcell withmultipleelectro-activelayersinthispaper. Regardingourpreviouswork,thismodelingmethod hasalreadyhelpeduswithimplementationofasensitivityanalysiswork[16]. Thismodelisappliedto revealthethermal-electric-electrochemicalinteractionsandtoillustratelocalheterogeneousevolutions ofphysicalprocessesindifferentlayersundervariousdischargeC-rateswithinthepouchcell.Toenable computationalefficientsimulation,asimplifiedelectrochemicalelectrodemodelonelectrodelevel isintroducedandcoupledtoa3Dthermal-electricmodeloncelllevel. Meanwhile, thesimplified electrochemicalmodelhasbeensuccessfullyvalidateduptoa4Cdischargeratebycomparingittoa full-orderpseudo-2D(P2D)model[17]. Inthiswork,wewillfullyintroduceourgoverningequations and geometries of the 3D multiphysics model at first, then the computational framework will be presented and discussed in detail. After a numerical validation of our electrochemical submodel, thethermal-electric-electrochemicalsimulationresultswillbedisplayedanddiscussedwithaspecial focusonheterogeneity. 2. 3DMultiphysicsModel The 3D multiphysics model uses a hierarchical framework to implement the coupling of the differentphysico-chemicalprocesseswithinthegeometryofthepouchcellpresentedinthiswork. Theframeworkandthecellgeometryforthis3DmodelareschematicallypresentedinFigure1. Theframeworkconsistsoftwosubmodelsontwolevels. Thesearetheelectrodeandthecelllevel. Oneachlevel,anindependentcoordinatesystemisassignedforthediscretization. Ontheelectrode level,thiscoordinatedescribesthelocalelectrochemicalprocessesofthecellalongthex-axisofeach electrochemicalelement,whichconsistsofanode,cathode,separator,andcurrentcollectors,intotalfive computationaldomains. Onthecelllevel,thecellconfigurationincludescurrentcollectortabs,current collectors,polymercover,andelectro-activelayersascomputationaldomainsinthreedimensions. Eachmeshedelementinanelectro-activelayerincludesallcomputationaldomainsoftheelectrode Energies2018,11,2998 3of26 levelandisdefinedasanonlinearresistor,whichrepresentsthelocalpotentialdropthatstemsfrom thelocalelectrochemicalprocesses. Onthecelllevel,chargeandenergyconservationequationsare solvedinthreedimensions. Bothsubmodelsarecoupledbyinter-levelvariables. Thesubmodelsand theircouplingapproachisprovidedindetailinthefollowingsections. Cell level Current collector tabs Current collector of anode (Cu) Current collector of cathode (Al) Electro-active layer Φ Φ Φ Φ Φ - + - + - ...... Polymer cover Cell ... Electro-active layers Y level X Z e- -Ilayer Ilayer Discharge Load Charge e- Separator Anode Cathode Li+ Li+ Y X 0 δan δsep δca x Z Electrode level Figure1.Illustrationofthehierarchicalframeworkin3DmultiphysicsmodelofaLithium-ionpouch cell,includinginformationofcellcomponents,computationaldomains,andcellgeometry. 2.1. ElectrodeLevel: ElectrochemicalSubmodel(ESM) In the electrochemical submodel (ESM), the electrochemical processes are modeled on the electrode level as shown in Figure 1, based on P2D models [17–19]. The discretization in both x andparticleradiusr directionsleadtohighcomputationalcost. Toreducethepresentcomplexity, twosimplificationscomparedtoafull-orderP2Dmodelareintroduced. InESM,thepore-wallfluxes j attheelectrode-electrolyteinterfacesi = {cathode,anode}aredenotedasaveragelumpedvariables x,i overspace[20,21]: 1 (cid:90) j¯ = j dx, i = ca,an (1) x,i x,i δi domain:i Energies2018,11,2998 4of26 anditcanbecalculatedaccordingto i j¯ = x , (2) x,i a Fδ i i with δ electrode thickness, i the local current density on the current collector flowing out of the i x electrode, F the Faraday constant and a the specific area of the electrode. The average lumped i variablej¯ simplifiesthelithiumdiffusionprocessintheactiveparticles,whichisusuallyexpressed x,i byFick’slaw: (cid:18) (cid:19) ∂c 1 ∂ ∂c s,i = r2D s,i , (3) ∂t r2 ∂r s,i ∂r withtheboundaryconditionsofthesurfacefluxattheouterboundary(r = R)ofparticlesandinthe i particlecenter(r =0)beinggivenasfollows: ∂c Ds,i ∂rs,i|r=Ri = j¯x,i (4) ∂c Ds,i ∂rs,i|r=0 =0 (5) whereR istheparticleradiusandD isthesoliddiffusioncoefficientoftheelectrode. Bytakinginto i s,i accounttheaforementionedboundaryconditions,Equation(3)canbeapproximatelysolvedasan analyticalsolutionofsurfaceandaveragelithiumconcentrations[22]. Thesurfaceconcentrationof lithiumineachelectrodeisexpressedas: cs,surf,i=cs,init,i−3jR¯x,iit− 2jD¯x,si,Rii (cid:32)110−n∑N=1λ12n(cid:33)1−e−λN+R12iDs,it+(cid:115)DπRs,i2iterfc(cid:32)λN+1(cid:115)DRs2i,it(cid:33) (6) − ∑N 2j¯x,iRi 1−e−λ2NR2iDs,it n=1λNDs,i wherec istheinitiallithiumconcentrationintheelectrodeparticle,Nisthetruncationerrorterm s,init,i number,andisdefinedasN =6inthiswork,andλ isthentheigenvalue,whichcanbecalculatedas n follows[22] λ −tan(λ ) =0, n =1,2,3,... (7) n n Thebulkconcentrationinthesolidphasesc¯ canbeevaluatedbysolvinglumpedmassbalance: s,i dc¯ j¯ s,i = −3 x,i, (8) dt R i Itisnotedthatthestateofcharge(SOC)ofthebatteryisalsorelatedtothebulksolidconcentration. Itisdefinedbythevariationrangeofbulksolidconcentrationincathodes[21],as c¯ −c¯ SOC= s,i s,init,i . (9) c −c s,0%,i s,init,i The initial solid concentration c is set as the value when SOC = 100% and the bulk solid s,init,i concentrationatSOC=0%,c canbeevaluatedby s,0%,i 12Ah c = c + . (10) s,0%,i s,init,i N ε Fδ A layer s,ca ca Energies2018,11,2998 5of26 Thesoliddiffusionoverpotentialineachelectrodeη isderivedasafunctionofsurfaceand diff,i bulkconcentrations,andisexpressedas: η = E(cs,surf,i)−E( c¯s,i ), (11) diff,i i c i c s,max,i s,max,i where c is the maximum lithium concentration in active particles, and E is the open circuit s,max,i i potential(OCP)ineachelectrode. TheButler-Volmerequationdescribesthechargetransferreaction rateattheelectrode-electrolyteinterfaces. Thereactionrateisdrivenbythekineticoverpotentialη , et,i whichisdefinedasthedeviationofthepotentialinsolidphaseΦ ,aswellasthepotentialinthe s,i solutionphaseΦ ,andE, e i (cid:20) (cid:18) (cid:19) (cid:18) (cid:19)(cid:21) i α Fη α Fη j¯ = 0,i exp a,i et,i −exp − c,i et,i , (12) x,i a F RT RT i c η = Φ −Φ −E( s,surf,i) (13) et,i s,i e i c s,max,i whereα andα arethechargetransfercoefficients,Ristheidealgasconstant,andi istheexchange a,i c,i 0,i currentdensity,whichisgivenasafunctionofthesurfacelithiumconcentrationinsolidphasesand lithiumconcentrationinthesolutionphaseattheelectrolyte-electrodeinterface. (Fortheevaluationof voltageinECMelement,allfieldvariables,e.g.,Φ ,c ,Φ ,intheanodeareevaluatedatx =0andin e e s thecathodeatx = L[21]). i = k Fcαa(cid:16)c −c (cid:17)αacαc (14) 0,i 0,i e s,max,i s,surf,i s,surf,i Here k is a kinetic rate coefficient, c is the lithium concentration in solution phase at the 0,i e electrolyte-electrodeinterfaces. Assumingthatα = α = α =0.5,basedonEquation(12),thekinetic a,i c,i overpotentialη canbeanalyticallysolved[21]: et,i RT (cid:18) (cid:113) (cid:19) η = ln Ψ + Ψ2+1 (15) et,i αF i i wherethetermΨ isgivenbyEquation(16): i Ri F Ψ = i x , (16) i 6ε i δ i 0,an i andε istheporosityofelectrodei. i Forsimplification,thevoltageineachESMelementcanbedenotedasthesumofdifferentterms, asdescribedbyPradaetal. [21]: (cid:18) (cid:19) (cid:18) (cid:19) c¯ c¯ Φ (L)−Φ (0) = E s,ca −E s,an +η −η +Φ (L)−Φ (0), (17) s,ca s,an ca an ca an e e c c s,max,ca s,max,an ThereferencepotentialontheanodeΦ (0)isdefinedbytheelectricpotentialatthenegative s,an currentcollectoroncelllevelΦ−. Thetotaloverpotentialηi isdefinedasthesumofoverpotentialsof diffusionandkinetics: η = η +η (18) i diff,i et,i TheelectricpotentialinelectrolyteΦ isevaluatedinthecontinuousliquidphase,whichconsists e of3domains,i = cathode,anode,separatorattheelectrodelevel: ∇·(cid:16)κeff∇Φ (cid:17)−∇(cid:34)κieffRT (cid:18)1+ dln f (cid:19)(cid:16)1−t0 (cid:17)∇lnc (cid:35)+a Fj¯ =0, (19) i e F dlnc + e i x,i e Energies2018,11,2998 6of26 wherej¯ =0. Theboundaryconditionsaregivenasfollows: x,sep ∇Φe|x=0,L =0 (20) andattheinterfaces, −−κκsaeeenffpff∇∇ΦΦee||xx==δ(a−δnan=+δs−epκ)−seefp=f∇−Φκece|axf=f∇δa+nΦe|x=(δan+δsep)+ (21) eff whereκ istheeffectiveionicconductivityofelectrolyteineverydomain,c issaltconcentrationin i e electrolyteandt0 isthetransferencenumberofLi+. WithEquations(19)–(21),Φ hasbeencalculated + e analyticallyineverydomainduetothelumpedreactionrateasfollows[21]: 2RT (cid:18) dln f (cid:19)(cid:16) (cid:17) c (x) i x2 Φ (x) = Φ (x =0)+ 1+ 1−t0 ln e − x (22) e e F dlnc + c (x =0) eff e e 2δ κ an an fortheanodedomain, 2RT (cid:18) dln f (cid:19)(cid:16) (cid:17) c (x) i δ i (x−δ ) Φ (x) = Φ (x =0)+ 1+ 1−t0 ln e − x an − x an (23) e e F dlnc + c (x =0) eff eff e e 2κ κ an sep fortheseparatordomain,and 2RT (cid:18) dln f (cid:19)(cid:16) (cid:17) c (x) i (L−x)2 Φ (x) = Φ (x =0)+ 1+ 1−t0 ln e + x e e F dlnc + c (x =0) eff e e 2δ κ ca ca (cid:32) (cid:33) (24) −ix δan +2δsep + δca 2 eff eff eff κ κ κ an sep ca forthecathodedomain. Theeffectiveionicconductivitiesaredefinedbyκeff =κεbruggi,andε isthe i i i porosityineverydomainandbrugg istheBruggemanncoefficient. i Themassconservationoflithiumsaltintheelectrolyteisgivenby: ε ∂ce = ∇·(cid:16)Deff∇c (cid:17)+(cid:16)1−t0 (cid:17)a j¯ , i = ca,an,sep (25) i ∂t e,i e + i x,i withthecorrespondingboundaryconditionsexpressedby: Dee,fif∇ce|x=0,L =0 (26) −−DDeeee,,ffsaeffnp∇∇ccee||xx==δ(a−nδan=+δ−sepD)−ee,fs=efp∇−cDe|eex,fc=afδ∇a+nce|x=(δan+δsep)+ (27) eff Here, D is the effective diffusion coefficient of electrolyte, which is evaluated by e,i Deff = D εbruggi for every domain. The governing equations of the electrochemical submodel, e,i e i all summarized parameters and material properties of the ESM are listed in Appendixes A and B,respectively. Furtherequationsabouttemperaturedependenceandthermalbehavioraregivenin thefollowingsections. 2.2. Thermal-ElectricSubmodel Onthecelllevel,a3Dthermal-electriccontinuumsubmodelisdevelopedtosolvetemperature T, and local electric potentials Φ+ and Φ− at the current collectors and tabs. A 3D geometry of a Energies2018,11,2998 7of26 lithium-ionpouchcellwithmultipleelectro-activelayersisconsideredtobeillustratedinFigure1. Itconsistsof6computationaldomains—twoelectrictabs,twoelectriccurrentcollectors,electro-active layers,andpolymercover.Thereare40parallellayersintheelectro-activedomain.Theheatconvection boundaryconditionsaresetonthedomainsofpolymercoverandtabs.Thermalandelectricproperties ofcellcomponentsareconsideredtobeanisotropicintheXY-planeandinZ-directions[7]. Theenergy conservationforalldomainsyields: ∂T ρC = ∇·(λ∇T)+q˙ , (28) p X ∂t wherethematerialdensityofcellcomponentρ,specificheatcapacityC andthermalconductivityλ p areattributedtoeverydomain. Thevolumetricheatsourceintheelectro-activedomainq˙ isthesum X ofthereactionheatq˙X,r,theentropychangeq˙X,S andJouleheatq˙X,Ω: q˙X = q˙X,r+q˙X,Ω+q˙X,S. (29) All locally volumetric heat sources are determined by the electrochemical submodel on the electrode level. The coupling values and expressions will be introduced in the following section. In present simulation scenarios, the cell surfaces are assumed to be exposed to the environment, wheretheconvectiveheattransferboundaryconditionsonsurfacesincludingtabsyield: −(cid:126)n·λ∇T = α (T−T ), (30) h amb where(cid:126)nisthenormalunitvector, T isconstantambienttemperatureandα istheheattransfer amb h coefficientbetweencellandtheenvironment. ThechargeconservationofthecellfollowsOhm’slaw. Intabandcurrentcollectordomains,itisexpressedas: ∇·(cid:0)σ∇Φ (cid:1)+i =0, j = −,+ (31) j j X whereσ istheelectricalconductivityofelectrodecurrentcollector(+or−),andi isthevolumetric j X currentdensityflowingfromelectro-activelayers. Intabdomains,i iszero. Foradischargeprocess, X theelectricboundaryconditionsasshowninFigure1areexpressedas I (cid:126)n·(−σ−∇Φ−) = (32) A− atthenegativeelectrodetab,and I (cid:126)n·(−σ+∇Φ+) = − (33) A+ atthepositiveelectrodetab,where A− and A+ areareasofanodeandcathodetabs,respectively. On thetabterminalofanodecurrentcollectors,thereferencepotentialisdefinedas: Φ− =0 (34) A1Dnonlinearresistornetworkisemployedforchargeconservationinelectro-activedomain insteadofPoisson’sequationinZ-direction[23]. Thenumberofnonlinearresistors N isequalto e thenumberofnodesontheXY-planeaftermeshingthecell(seeFigure1). AccordingtoGerver’s approach[23],thelocalcurrent I intheelectro-activedomainfollowstherelation: node Φ −Φ I = +,node −,node (35) node R node whereΦ+,node andΦ−,node arethelocalpotentialsatnodesoncathodeandanodecurrentcollectors ontheXY-plane,respectively. R (I ,η )isthelocalresistanceevaluatedbyanindependent node node node Energies2018,11,2998 8of26 electrochemicalsubmodel. Usingnonlinearresistors,itispossibletocombinelocalelectrochemical processes withgeometryonthe celllevel. Therefore, thecouplingof toomanypartialdifferential equations(PDEs)ofelectrochemicalprocessesisavoided,andthesimulationofchargeconservationin 3Dpouchcellconfigurationwithafasterconvergencethanfull-physicscouplingisenabled[23]. 2.3. ModelCoupling Figure1showsthecouplingofthetwosubmodels. Eachindependentelectrochemicalsubmodel performsasanonlinearresistoroncelllevel. Inthefollowingwewillintroducethecouplingstrategy betweenbothlevelsubmodelsandthecomputationalframework. 2.3.1. Inter-LevelCouplingVariables In this 3D multiphysics model, variables on the electrode level are firstly solved, then these solutionsaredeliveredtothecelllevelthroughinter-levelcouplingvariablesasshowninFigure2. Intheelectrochemicalmodel,localresistanceR (X,Y,Z)andvolumetricheatsourceq˙ (X,Y,Z)are node X determined. Thethermal-electricsubmodelprovideslocaltemperatureT(X,Y,Z),electricpotentials Φ+(X,Y,Z) and Φ−(X,Y,Z), and local current Inode(X,Y,Z) as an input for the electrochemical submodels. Thelocaltemperatureofanelectro-activeelementistheaveragevalueofitscorresponding nodes. ThelocalresistanceR isdefinedasafunctionoflocalcurrent I andlocalpotentialdrop node node betweenanodeandcathodecurrentcollectors. Thepotentialdropiscoupledonbothlevels: Φ+,node−Φ−,node = Φs,ca(L)−Φs,an(0) (36) ThenusingEquations(15)–(17),(22)and(24),thispotentialdropcanbefinallywrittenas: Φ+,node−Φ−,node =Eca(θca)−Ean(θan) (cid:32) (cid:112) (cid:33) RT Ψ + Ψ2 +1 + ln ca (cid:112) ca αF Ψ + Ψ2 +1 an an 2RT (cid:18) dln f (cid:19)(cid:16) (cid:17) c (L) (37) + 1+ 1−t0 ln e F dlnc + c (0) e e (cid:32) (cid:33) − ix δan +2δsep + δca . 2 eff eff eff κ κ κ an sep ca whereθ = c /c . I isalocalcurrentofthecorrespondingresistorandisdenotedby: i s,surf,i s,max,i node i A I = x (38) node N elem whereN isthetotalnumberofelectrochemicalsubmodelelementsinanelectro-activelayer. elem Usingasimilarapproachforthecouplingofelectrochemicalprocessesandchargeconservationon bothlevels,thevolumetricheatsourcesq˙X,r,q˙X,Ω andq˙X,S aredefinedasaveragevaluesonelectrode levelasfollows[24]: q˙ = ∑ a Fj¯ η, (39) X,r i x,i i i=ca,an q˙X,Ω = aiFj¯x,i[Φe(L)−Φe(0)] (40) and q˙ = ∑ a Fj¯ T∂Ei (41) X,S i x,i ∂T i=ca,an Usingthethreetermsabove,thevolumetricheatsourceq˙ ofeverycorrespondingnodeinthe X resistorelementcanbeevaluatedbyEquation(29). Forexample, iftheresistorisalinearelement Energies2018,11,2998 9of26 withtwonodes,thenonenoderepresentstheanode,andtheotheronedenotesthecathode. Also, everyvolumetricheatsourceisdividedintotwopieces:anodeandcathode,andthentheyareassigned accordingly. For the ohmic heat loss, it is assumed that the active material is well electronically conductive,thustheohmiclossinsolidphasesisneglectedhere. Electrode level Cell level . R q R (t) node X Anode Cathode . node q (t) X c (t) s,Ri η(t) i Li+ Y Φ (x,t) s,i Φ(x,t) e X Z c(x,t) e T Φ Φ I - + node T(X,Y,Z,t) I (X,Y,Z,t) node Separator 0 x Φ-(X,Y,Z,t) Φ+(X,Y,Z,t) Figure2.Overviewoftheinter-levelcouplingvariablesbetweentheelectrodelevelandthecelllevel inthe3Dmultiphysicsmodel. 2.3.2. ComputationalFramework ReferringtoAllu’s[14]andKim’s[15]works,theframeworkofthe3Dmultiphysicsmodelwith theinter-levelcouplingvariablesenablesdifferentsolverstobedeployedwiththedifferentsubmodels. Figure3illustratesthecouplingalgorithmindetail. Picard iteration Recent state @ t Refresh t = t i i+1 i Element #1 … Electrochemical submodel Element #N Δt Self-consistant iteration T, Φ, Φ , I - + node 3D thermal-electric Electrochemical state @t i+1 state @t i+1 . q R X node Δt 3D thermal-electric Self-consistant submodel iteration Picard iteration Refresh t = t Recent state @ t i+1 i i Figure3.Illustrationofthecomputationalalgorithmforthe3Dmultiphysicsmodel. Energies2018,11,2998 10of26 The3Dmultiphysicsdiagramshowsthattheelectrochemicalsubmodelprovidesandcalculates nodalandelementalelectrochemicalinformation q˙ and R intheelectro-activedomainonthe X node celllevel,andthe3Dcell-levelsubmodelprovidesaveragedlumpedinputs T, Φ and I forthe j node electrode-levelmodel. Atthebeginningofeachsimulation,onbothelectrodeandcelllevels,initial conditions are evaluated by static calculations of the two submodels. After starting the transient simulation,atthestartofanytimestep,t,allstoredsolutionsandinter-levelcouplingvariableson i bothlevelsaredefinedbytherecentstateatt. Themeshingresolutionandthenumberofnodesinthe i electro-activedomaindeterminethetotalnumberofelectrochemicalsubmodels. Allelectrochemical submodelsreadtheirinter-levelcouplingvariablesfromthe3Dthermal-electricsubmodelandcalculate the followingstates at ti+1 with aself-consistent iteration. Consequently, the exchangestate of all correspondingelectrochemicalsubmodelsat ti+1 andtherecentstateofthe3Dsubmodelat ti are employedinthe3Dthermal-electricsubmodeltocalculateitsnewexchangestatesatti+1. Attheend ofthisiterativeprocess, newexchangestatesofbothsubmodelsbecometherecentstates, andthe 3Dmultiphysicsmodelcontinueswiththenextglobaltimestep. Thisiterativeapproachofthe3D multiphysicsmodelisbasedonPicarditerativemethod,whichhasbeenusedfortransientcalculation ofmultiphysicsmodels. Itenableseachsubmodeltoproceedwithitsownself-consistentiterationsto reachaconvergencecriterionaftereachglobaltimestepisincreased[14]. 2.4. ModelStudies In this paper, the electrochemical submodel is validated first at various current rates with a full-orderP2Dmodel[17],thenthe3Dmultiphysicsmodelisusedtosimulatea12Ahlarge-format pouchcellatvariousdischargeratesundergalvanostaticdischargeconditions. Theelectrodesofthe cellarecomposedoflithiumrichNCMmaterials(Ni,Co,Mn)andgraphite,respectively. Weassumea binaryelectrolytecontainingLiPF inethylenecarbonate/ethylmethylcarbonate. Thecellproperties 6 are presented in the Appendixes B and C as model parameters. The cell dimensions were 99 mm width × 120mm height × 9 mm thickness. The cell contains 40 parallel stacked electro-active layers, regarding the normal number of layers for commercial pouch cell [12]. The multiphysics modelwasdefinedasanopensystem,withaconvectiveheatexchangecoefficientof15W/m2·K connected to the environment. The initial and ambient temperatures were chosen as 25 ◦C. The electrochemicalsubmodelswereimplementedinMATLAB(R2015b,MathWorks,Natick,MA,USA) with SundialsTB solvers [25]. The computation cost of electrochemical submodels is estimated by MATLABusingthefunctiontic-toc;itcanrecordtheinternalexecutiontimeforMATLABprograms. The 3D thermal-electric submodel was implemented in Ansys APDL 15.0. All simulations were executedonan8coreI7-2600processorwith16GBmemory. 3. ResultsandDiscussion 3.1. NumericalValidationofElectrochemicalSubmodel TheESMisvalidatedbycomparingthedischargebehaviorformultipleC-ratestoafull-order P2Dmodel[17]asshowninFigure4. Figure4aillustratesthatthedischargecurveofESMisingoodagreementwiththefull-order P2Dmodelbelow4Cdischarge,butitshowsalargerdeviationat4C.TherelativeerrorofESMat 4Cfordischargevoltagenearlyincreasesto5%intheendasshowninFigure4b. Forthesimulation carriedoutinthiswork,thiserrorisstillwithinanacceptablerange,becauseitdoesnotleadtolarge deviationsofinter-levelcouplingvariablesforthemultiphysicsmodel. Thedeviationislikelycaused bytheassumptionoflumpedpore-wallflux,whichresultsinoverestimatedreactionratesalongthe anode. Asaconsequence,ESMoverestimatestheconsumptionofLi+ attheanodeandleadingto higher overpotentials. Nevertheless, the final capacity is identical, which suggests that the entire processisnotaffectedbythisoverestimation.