Multiphysics Modeling for Detailed Analysis of Multi-Layer Lithium-Ion Pouch Cells

26 Pages·2017·1.37 MB·English

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.:[phone]
(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.