基于統(tǒng)一相場(chǎng)理論的鋰電池電極顆粒斷裂模擬研究1)

吳建營(yíng) 洪屹峰

* (華南理工大學(xué)亞熱帶建筑科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,廣州 510641)

? (華南理工大學(xué)土木工程系,廣州 510641)

引言

隨著化石能源的日益枯竭,清潔可再生的新能源開(kāi)發(fā)和穩(wěn)定可靠的能量?jī)?chǔ)存受到人們的密切關(guān)注.其中,鋰離子電池由于其高容量、長(zhǎng)循環(huán)壽命等優(yōu)點(diǎn)被廣泛應(yīng)用于電子產(chǎn)品、電動(dòng)汽車(chē)、電網(wǎng)等領(lǐng)域[1].由此,也對(duì)鋰離子電池的失效分析和控制管理提出了更高的要求[2].在眾多導(dǎo)致鋰離子失效的因素中,力學(xué)行為的影響至關(guān)重要[3-4],尤其是,在充放電過(guò)程中電極顆粒表面發(fā)生鋰離子的脫出或嵌入,電極顆粒材料內(nèi)部發(fā)生不均勻的體積變化,從而產(chǎn)生拉應(yīng)力并導(dǎo)致電極材料開(kāi)裂和破壞[5-6].研究表明,電極顆粒的力學(xué)失效將會(huì)對(duì)鋰離子電池的性能產(chǎn)生嚴(yán)重不利影響,具體包括: 電子電導(dǎo)率衰減[7]、電池容量降低[8]以及循環(huán)壽命折減[9]等.因此,研究并預(yù)測(cè)鋰電池電極顆粒在充放電過(guò)程中的力學(xué)失效行為,對(duì)開(kāi)發(fā)具有高能量密度、高安全性能和長(zhǎng)循環(huán)壽命的鋰電池具有重要意義.

近年來(lái),國(guó)內(nèi)外學(xué)者針對(duì)鋰離子電池電極顆粒在充放電過(guò)程中的力學(xué)失效行為展開(kāi)了廣泛的實(shí)驗(yàn)和理論研究.Fincher 等[10]測(cè)試了金屬鋰負(fù)極材料的力學(xué)性能;Zhang 等[11]基于變分原理構(gòu)建了擴(kuò)散–力學(xué)耦合塑性–損傷模型,研究了將硅負(fù)極材料的力學(xué)行為.Xu 等[12]通過(guò)試驗(yàn)研究了無(wú)機(jī)鋰鹽化合物陽(yáng)極材料在充放電循環(huán)過(guò)程中的力學(xué)失效行為；Klinsmann 等[13-14]等學(xué)者采用脆性斷裂相場(chǎng)模型,研究了嵌鋰和脫鋰過(guò)程中化學(xué)擴(kuò)散–力學(xué)變形耦合引起的陽(yáng)極顆粒開(kāi)裂行為;Zhang 等[15]則采用內(nèi)聚裂縫模型,模擬了在充放電循環(huán)過(guò)程下電極顆粒的裂縫發(fā)展.

上述工作為研究鋰電池電極顆粒的力學(xué)失效行為奠定了良好的基礎(chǔ).然而,遺憾的是,這些研究大多采用內(nèi)聚裂縫模型或脆性斷裂相場(chǎng)模型: 前者假定裂縫尖端的最大應(yīng)力為一有限大小的材料參數(shù)(即內(nèi)聚力強(qiáng)度或破壞強(qiáng)度),并采用基于強(qiáng)度的起裂準(zhǔn)則描述裂縫起裂,雖然可以引入裂縫方向判據(jù)以確定裂縫擴(kuò)展路徑,但仍然難以考慮多裂縫、任意擴(kuò)展路徑等復(fù)雜行為,且數(shù)值模擬結(jié)果存在難以解決的網(wǎng)格敏感性問(wèn)題;對(duì)于后者,無(wú)需如Griffith 線彈性力學(xué)般預(yù)設(shè)初始裂縫,且自帶基于能量的裂縫擴(kuò)展判據(jù)和基于系統(tǒng)能量最小的裂縫方向判據(jù),代表了線彈性斷裂力學(xué)的最新成果,然而其分析結(jié)果存在嚴(yán)重的裂縫尺度敏感性問(wèn)題: 裂縫尺度越小,預(yù)測(cè)給出的峰值應(yīng)力越高,因此無(wú)法完全解決裂縫起裂問(wèn)題[16-17],特別是難以描述完好或僅有弱奇異性固體的裂縫演化特征.

2017 年以來(lái),Wu[18]將斷裂力學(xué)與損傷力學(xué)相結(jié)合,建立了同時(shí)適用于脆性斷裂和準(zhǔn)脆性破壞的統(tǒng)一相場(chǎng)理論.從該理論出發(fā),不僅脆性斷裂相場(chǎng)模型[19-20]可作為其特例直接給出,還首次提出了一類(lèi)相場(chǎng)內(nèi)聚裂縫模型即PF-CZM[21-22].作為一種裂縫正則化變分方法,相場(chǎng)內(nèi)聚裂縫模型將基于強(qiáng)度的裂縫起裂準(zhǔn)則、基于能量的裂縫擴(kuò)展準(zhǔn)則以及基于變分原理的裂縫路徑判據(jù)納入同一框架內(nèi).該模型非常便于通過(guò)有限元等方法加以數(shù)值實(shí)現(xiàn),大量算例均表明其分析結(jié)果不依賴(lài)于裂縫尺度和有限元網(wǎng)格等數(shù)值參數(shù),且僅需少量標(biāo)準(zhǔn)材料參數(shù)即可準(zhǔn)確描述固體結(jié)構(gòu)的損傷與破壞行為,具有十分優(yōu)異的預(yù)測(cè)能力.歸因于上述優(yōu)點(diǎn),統(tǒng)一相場(chǎng)理論和相場(chǎng)內(nèi)聚裂縫模型一經(jīng)提出,就迅速得到了國(guó)內(nèi)外學(xué)者的廣泛認(rèn)可和直接應(yīng)用[23–27].

為此,本工作擬在統(tǒng)一相場(chǎng)理論框架內(nèi),建立考慮化學(xué)擴(kuò)散–力學(xué)變形耦合影響的相場(chǎng)內(nèi)聚裂縫模型,發(fā)展相應(yīng)地多場(chǎng)有限元數(shù)值實(shí)現(xiàn)算法,并應(yīng)用于鋰電池電極顆粒的力學(xué)失效行為研究.一方面,通過(guò)復(fù)現(xiàn)一些前人相關(guān)研究中的結(jié)果,驗(yàn)證統(tǒng)一相場(chǎng)理論擴(kuò)展應(yīng)用于鋰電池電極顆粒開(kāi)裂引發(fā)力學(xué)失效問(wèn)題的適用性.另一方面,克服此前基于脆性斷裂相場(chǎng)模型相關(guān)研究?jī)H考慮帶有初始缺陷電極顆粒的缺點(diǎn),預(yù)測(cè)完好電極顆粒裂縫演化導(dǎo)致的力學(xué)失效.在此基礎(chǔ)上,以期進(jìn)一步發(fā)現(xiàn)更多潛在裂縫模式,為鋰電池電極顆粒的優(yōu)化設(shè)計(jì)提供科學(xué)依據(jù).

1 化學(xué)–力學(xué)耦合相場(chǎng)內(nèi)聚裂縫模型

不失一般性,考慮內(nèi)嵌裂縫或界面的固體Ω ?Rndim(ndim=1,2,3),其外邊界和外法向矢量分別記為 ?Ω 和n,材料點(diǎn)的空間坐標(biāo)記為x.對(duì)于力學(xué)場(chǎng),邊界 ?Ω 劃分為互補(bǔ)的兩部分 ?Ωu和 ?Ωt,分別施加給定的位移邊界u*(x) 和力邊界t*(x) ;對(duì)于化學(xué)場(chǎng),則將邊界 ?Ω 劃分為 ?ΩJ和 ?ΩC兩部分,分別施加離子通量邊界J*(x) 和離子濃度邊界C*(x).

雖然統(tǒng)一相場(chǎng)理論也適用于有限變形[28],這里僅考慮小應(yīng)變情況,故采用位移場(chǎng)u(x) 和線性應(yīng)變場(chǎng) ε(x)=?su(x) 描述固體的變形狀態(tài),這里 ?s(·) 表示對(duì)稱(chēng)梯度算子.伴隨鋰電池充放電過(guò)程中電極顆粒表面鋰離子的脫出或嵌入,離子濃度C(x) 發(fā)生改變,導(dǎo)致電極顆粒材料內(nèi)部呈現(xiàn)不均勻的體積變化,通常采用附加化學(xué)應(yīng)變 εc加以描述,并假定其大小與離子濃度變化成正比.由此,將應(yīng)變張量 ε 分解為彈性應(yīng)變張量 εe和化學(xué)應(yīng)變張量 εc之和的形式,即

式中,VH代表鋰離子的偏摩爾體積,C0為鋰離子的初始濃度;I為二階單位張量.

相場(chǎng)裂縫模型將尖銳裂縫/界面 S 彌散為裂縫帶B ?Ω,其尺度記為b＞0,并引入裂縫相場(chǎng)d(x)∈[0,1]及其空間梯度 ?d(x) 描述其狀態(tài).裂縫帶 B 的外邊界記為 ?B,外法向矢量表示為nB.需要指出,這里裂縫帶 B 并非預(yù)設(shè)或固定,而是遵循自身特定的本構(gòu)關(guān)系發(fā)生起裂、擴(kuò)展和演化等過(guò)程.類(lèi)似地,裂縫相場(chǎng)也可以施加某種強(qiáng)制邊界條件.例如,對(duì)于彈性區(qū)域有d(x)=0,對(duì)于初始裂縫有d(x)=1 等.

1.1 擴(kuò)散–變形耦合過(guò)程的控制方程

對(duì)于電極顆粒表面鋰離子脫出或嵌入引起的化學(xué)擴(kuò)散–力學(xué)變形過(guò)程,其虛功原理表示為

式中,b*為均布體積力,ρ 和分別為材料密度和位移加速度;J為離子通量,μ 為描述離子擴(kuò)散的化學(xué)勢(shì).利用高斯–斯托克斯散度定理,可以給出

上述離子濃度和位移子問(wèn)題相互耦合,共同構(gòu)成化學(xué)擴(kuò)散–力學(xué)變形過(guò)程的控制方程.

同時(shí),若忽略上述耦合過(guò)程中的溫度變化(即恒溫假定)并考慮控制方程(3),熱力學(xué)第二定律表示為如下Clausius–Duhem 不等式

這里 ψ(ε,C,d) 為固體的Helmholtz 自由能勢(shì).由此,可以給出如下一般形式的化學(xué)和力學(xué)本構(gòu)關(guān)系

相應(yīng)地,能量耗散不等式退化為

式中,化學(xué)耗散 D˙c和力學(xué)耗散 D˙m分別表示為

這里假設(shè)化學(xué)耗散和力學(xué)耗散相互解耦且均為非負(fù)值,以保證在純化學(xué)擴(kuò)散或純力學(xué)變形情況下能量耗散不等式均得到滿足.

離子通量J由修正的菲克定律給出[29-30]

這里,M(C) 為離子遷移率,D為擴(kuò)散系數(shù),C?=C/Cmax為歸一化離子濃度,Cmax為電極顆粒能容納的鋰離子最大濃度;R為理想氣體常數(shù),T為絕對(duì)溫度.可以看出,上述離子通量自動(dòng)滿足化學(xué)耗散不等式(7a).

1.2 裂縫相場(chǎng)演化的能量原理

根據(jù)圖像分割理論[31],尖銳裂縫的面積AS可以表示為如下體積分形式[32]

式中,裂縫面積密度函數(shù) γ(d;?d) 表示為[18]

這里,裂縫幾何函數(shù) α(d)∈[0,1] 是裂縫相場(chǎng)d(x) 的遞增函數(shù).可以嚴(yán)格證明: 當(dāng)相場(chǎng)尺度參數(shù)b趨近于零時(shí),式(9)給出的Ad(d) 收斂于尖銳裂縫的面積AS.

統(tǒng)一相場(chǎng)理論中,裂縫相場(chǎng)演化遵循如下能量準(zhǔn)則[18,21]

這里,材料斷裂能常數(shù)Gf表征單位面積裂縫擴(kuò)展所需的能量.從上式出發(fā)并結(jié)合高斯–斯托克斯散度定理,可以得到(具體推導(dǎo)過(guò)程詳見(jiàn)文獻(xiàn)[33]第1.2 節(jié))

式中,熱力學(xué)廣義力q和源項(xiàng)Q分別表示為

式(12)和式(13)構(gòu)成了裂縫相場(chǎng)演化的控制方程[18,33].

1.3 化學(xué)–力學(xué)耦合本構(gòu)關(guān)系

類(lèi)似于應(yīng)變分解式 (1),鋰電池充放電過(guò)程中電極顆粒材料的化學(xué)–力學(xué)耦合Helmholtz 自由能函數(shù)可以表示為

彈性勢(shì)能 ψe和化學(xué)勢(shì)能 ψc分別定義為

其中,E0為材料彈性剛度張量; ω(d)∈[0,1] 為單調(diào)遞減的能量退化函數(shù),表征裂縫演化對(duì)彈性應(yīng)變能的影響[18]; μ0為電極材料的初始參考化學(xué)勢(shì).需要說(shuō)明的是,式(15 a)假定彈性剛度張量E0與鋰離子濃度C無(wú)關(guān);如有必要,可以考慮鋰離子濃度C對(duì)彈性模量等材料力學(xué)參數(shù)的影響.

相應(yīng)地,由式(5)和式(8)給出化學(xué)擴(kuò)散和力學(xué)變形過(guò)程的本構(gòu)關(guān)系如下

上述力學(xué)本構(gòu)關(guān)系式(16b)和式(17)給出的材料受拉和受壓行為完全對(duì)稱(chēng),這與脆性斷裂和準(zhǔn)脆性破壞的特征不符.解決這一問(wèn)題最簡(jiǎn)單的策略是引入等效應(yīng)力將能量釋放率式(17b)重新定義為

這里,E0為材料彈性模量;是能量釋放率的初始閾值,與材料破壞強(qiáng)度f(wàn)t有關(guān);MacAuley 括號(hào)定義為 〈x〉=max(x,0).對(duì)于脆性斷裂和準(zhǔn)脆性破壞,等效應(yīng)力可取為有效應(yīng)力張量的最大特征值

另一方面,由于裂縫擴(kuò)展的不可逆特性,相場(chǎng)演化控制方程(12)為不等式.處理這一問(wèn)題最簡(jiǎn)單的方法是考慮加載歷史的影響,將式(13)中的裂縫相場(chǎng)源Q修改為

需要說(shuō)明的是,對(duì)于裂縫高速動(dòng)態(tài)擴(kuò)展問(wèn)題,可以通過(guò)哈密頓原理在裂縫相場(chǎng)控制方程(12a)中引入與裂縫相場(chǎng)一階時(shí)間導(dǎo)數(shù)(裂縫擴(kuò)展速率)有關(guān)的微阻尼或(和)與裂縫相場(chǎng)二階時(shí)間導(dǎo)數(shù)(裂縫擴(kuò)展加速度)有關(guān)的微慣性,以描述裂縫高速動(dòng)態(tài)擴(kuò)展時(shí)的能量耗散和應(yīng)變率效應(yīng),具體參見(jiàn)前期相關(guān)工作[34].對(duì)于鋰電池充放電過(guò)程脫嵌鋰引起的電極顆粒開(kāi)裂問(wèn)題,其應(yīng)變率效應(yīng)很小、也很難量化,因此這里忽略裂縫相場(chǎng)的時(shí)間導(dǎo)數(shù),僅考慮與位移速度有關(guān)的宏觀慣性.然而,這并非意味著鋰離子擴(kuò)散和裂縫擴(kuò)展相互解耦.事實(shí)上,相場(chǎng)控制方程中,裂縫驅(qū)動(dòng)力(為彈性應(yīng)變的函數(shù))與鋰離子濃度有關(guān),而后者則與時(shí)間有關(guān),自然考慮了鋰離子擴(kuò)散對(duì)裂縫擴(kuò)展的時(shí)間效應(yīng);反之,雖然裂縫擴(kuò)展在極短時(shí)間內(nèi)完成,但裂縫相場(chǎng)也會(huì)直接通過(guò)裂縫尖端的靜水應(yīng)力對(duì)鋰離子擴(kuò)散產(chǎn)生影響.因此,如圖1 所示,這里提出的模型考慮了離子擴(kuò)散過(guò)程和裂縫演化過(guò)程的雙向耦合.

圖1 鋰電池電極顆粒的擴(kuò)散-開(kāi)裂-變形耦合過(guò)程Fig.1 Chemo-mechanical coupling in storage particles of LIBs

1.4 脆性斷裂的本構(gòu)特征函數(shù)

已有研究表明[18,22]: 對(duì)于脆性破壞,統(tǒng)一相場(chǎng)理論可以采用如下裂縫幾何函數(shù)和能量退化函數(shù)

其中參數(shù)a1=4lch/(πb),這里Irwin 特征長(zhǎng)度lch:=表征材料脆性程度: 其值越小,材料越脆.

如圖2 所示,上述本構(gòu)特征函數(shù)將給出具有線性軟化曲線的相場(chǎng)內(nèi)聚裂縫模型[21],可用于描述脆性材料的損傷破壞過(guò)程.

圖2 具有線性軟化曲線的相場(chǎng)內(nèi)聚裂縫模型Fig.2 Phase-field cohesive zone model with a linear softening law

2 數(shù)值實(shí)現(xiàn)和求解算法

偏微分方程(3)和(12)構(gòu)成了如圖1 所示的擴(kuò)散–變形–開(kāi)裂過(guò)程的耦合控制方程,通常采用多場(chǎng)有限元方法加以求解.

2.1 控制方程的弱形式

為便于有限元數(shù)值實(shí)現(xiàn),將上述偏微分控制方程轉(zhuǎn)化為如下弱形式

式中,δu,δC和 δd分別為任意容許的節(jié)點(diǎn)位移、離子濃度和裂縫相場(chǎng)(變分);外力虛功和外邊界離子通量虛功分別表示為

需要說(shuō)明的是,這里采用式(19)給出的裂縫相場(chǎng)源項(xiàng),故控制方程(12)及其弱形式(21)為等式.

2.2 多場(chǎng)有限元數(shù)值實(shí)現(xiàn)

有限元空間離散時(shí),可將整個(gè)區(qū)域內(nèi)所有單元節(jié)點(diǎn)均賦予位移、濃度和裂縫相場(chǎng)自由度;為了減少計(jì)算量,也可將計(jì)算域 Ω 預(yù)先劃分為兩部分[35]: 可能出現(xiàn)裂縫的子區(qū)域 B 和剩余彈性部分 ΩB.對(duì)于前者,單元節(jié)點(diǎn)同時(shí)具有位移、濃度和裂縫相場(chǎng)自由度;而后者僅賦予位移和濃度自由度.已有研究表明[18,21-22,36-37],為保證計(jì)算精度,相場(chǎng)子區(qū)域 B 內(nèi)的單元大小一般取h≤b/5.

有限元空間離散后,位移、濃度和裂縫相場(chǎng)及其空間(對(duì)稱(chēng))梯度分別表示為(數(shù)值實(shí)現(xiàn)中,所有張量均采用Voigt 標(biāo)記)

于是,由弱形式 (21)可以給出如下殘量形式的平衡方程

需要指出,式(16 a)中的離子通量J涉及靜水應(yīng)力梯度 ?σH,在式(24b)數(shù)值實(shí)現(xiàn)時(shí)需要將單元積分點(diǎn)處的靜水應(yīng)力映射至單元節(jié)點(diǎn),具體詳見(jiàn)文獻(xiàn)[38].

2.3 時(shí)間域離散

方程組 (24)通常采用增量迭代法進(jìn)行求解,即將時(shí)間步 [0,t] 離散為N個(gè)增量步這里n=0,1,2,···,N-1.對(duì)于時(shí)間增量為Δt=tn+1-tn的典型子步在tn時(shí)刻所有狀態(tài)變量已知的條件下,求解tn+1時(shí)刻的狀態(tài)變量.

對(duì)于節(jié)點(diǎn)擴(kuò)散速率列向量a˙?,采用無(wú)條件穩(wěn)定的向后歐拉方法進(jìn)行時(shí)間域的離散,即

由此,將式 (24b)改寫(xiě)為

這里(下同)忽略了tn+1時(shí)刻物理量的下標(biāo)n+1.

節(jié)點(diǎn)位移a和速度a˙ 通常由Newmark 方法[39]給出

其中,節(jié)點(diǎn)試算位移atrial和試算速度a˙trial表示為

這里,系數(shù) β1∈[0,0.5] 和 β2∈[0,1] 對(duì)應(yīng)特定的時(shí)間積分方法.例如,β1=1/4 和β2=1/2 對(duì)應(yīng)隱式中點(diǎn)(梯形)方法; β1=0 和 β2=1/2 對(duì)應(yīng)顯式中心差分方法.

后續(xù)數(shù)值模擬均采用隱式Newmark 方法即 β1≠0.此時(shí),式(24a)可表示為節(jié)點(diǎn)位移an+1的殘量方程

2.4 交錯(cuò)迭代算法

式(24c)、式(27)和式(30)共同構(gòu)成了以節(jié)點(diǎn)裂縫相場(chǎng)、濃度和位移a為未知量的非線性代數(shù)方程組,通常采用迭代方法進(jìn)行求解.為了保證數(shù)值求解的健壯性,這里采用如圖3 所示的交錯(cuò)迭代算法.

圖3 求解濃度–相場(chǎng)–位移耦合控制方程的交錯(cuò)迭代算法流程圖Fig.3 Flow chart of the staggered algorithm in solving the coupled governing equations

(2)根據(jù)式(30)迭代求解節(jié)點(diǎn)位移a(k+1),此時(shí)節(jié)點(diǎn)濃度固定為更新后的結(jié)果而節(jié)點(diǎn)裂縫相場(chǎng)仍然取為上一迭代步給出的結(jié)果

(3) 基于更新后的節(jié)點(diǎn)位移a(k+1)和濃度由式(24c)迭代求解節(jié)點(diǎn)裂縫相場(chǎng)

當(dāng)控制方程的殘量范數(shù)小于某一給定容差,判斷為結(jié)果收斂;否則,繼續(xù)重復(fù)上述求解過(guò)程直至最終收斂.

已有數(shù)值計(jì)算經(jīng)驗(yàn)表明[18,21,33],上述交錯(cuò)迭代算法的收斂速度較慢,但數(shù)值健壯性極好,對(duì)于求解裂縫快速擴(kuò)展引起的脆性斷裂問(wèn)題十分有效.

3 電極顆粒力學(xué)失效分析應(yīng)用

上述化學(xué)-力學(xué)耦合相場(chǎng)內(nèi)聚裂縫模型適用于鋰電池充電(鋰離子脫出)、放電(鋰離子嵌入)過(guò)程中擴(kuò)散應(yīng)力引起的電極顆粒裂縫演化分析.由于篇幅受限,這里僅討論鋰離子脫出過(guò)程,鋰離子嵌入導(dǎo)致的電極顆粒力學(xué)失效分析將另文給出.

研究對(duì)象為錳酸鋰(LiMn2O4)鋰電池陽(yáng)極顆粒,材料參數(shù)如表1 所示[29,41-42,14].根據(jù)實(shí)際情況,考慮(圓或橢圓)柱體和球體兩種形狀的電極顆粒,分別采用二維和三維有限元進(jìn)行數(shù)值分析.與已有研究[13-14]一致,化學(xué)擴(kuò)散采用如下通量邊界條件

表1 錳酸鋰正極顆粒材料參數(shù)取值Table 1 Material parameters of the storage particle.

式中,c為表征鋰電池充放電快慢的參數(shù):c值越大,充放電速率越快;V和S分別表示電極顆粒的體積和表面積.鋰離子脫出時(shí),初始條件假定為C0=Cmax,邊界離子濃度達(dá)到最小值0 后保持恒定.

3.1 柱體電極顆粒: 二維分析

首先考慮如圖4 所示的圓形截面柱體電極顆粒,其截面半徑記為r0,此時(shí)有V/S=r0/2.

圖4 二維電極顆粒的幾何及邊界條件Fig.4 A cylindrical storage particle: Geometry and boundary conditions

為減少計(jì)算量,僅考慮圓柱體電極顆粒的截面行為,故可簡(jiǎn)化為二維問(wèn)題并假定為平面應(yīng)力狀態(tài).在某些情況下,還將根據(jù)問(wèn)題的對(duì)稱(chēng)性采用 1/2 或1/4模型進(jìn)行分析.由于模型計(jì)算結(jié)果與相場(chǎng)尺度參數(shù)b無(wú)關(guān),所有算例中均取為b=0.05 μm,網(wǎng)格劃分采用四節(jié)點(diǎn)雙線性插值單元.

充電過(guò)程中,隨著鋰離子逐漸脫出,電極顆粒表面濃度低于中心區(qū)域.此時(shí),式(1)中的附加化學(xué)應(yīng)變大小為負(fù)值,即電極顆粒表面發(fā)生收縮變形而產(chǎn)生環(huán)向拉應(yīng)力,而電極顆粒內(nèi)部產(chǎn)生環(huán)向壓應(yīng)力,如圖5 所示.

圖5 脫鋰過(guò)程電極顆粒的應(yīng)力狀態(tài)Fig.5 The mechanical stress state in the particle during Li extraction

3.1.1 彈性分析結(jié)果

首先考慮半徑為r0=5 μm 的電極顆粒,其充放電速率c=1,分析時(shí)間為t=5000 s.圖6 給出了最大主應(yīng)力隨時(shí)間的變化曲線.

圖6 電極顆粒彈性分析結(jié)果: 最大主應(yīng)力隨時(shí)間的變化曲線Fig.6 Elastic results of the particle: Evolution of the maximum principal stress

可以看出,上述電極顆粒的最大主應(yīng)力不超過(guò)0.30 MPa,明顯小于材料破壞強(qiáng)度f(wàn)t=1.3 GPa,此時(shí)材料處于線彈性狀態(tài).

在線彈性階段,附錄給出的解析解表明,圓柱體電極顆粒靜水應(yīng)力的拉壓分界線 σH(r)=0 與截面平均濃度等值線C(r)=重合.為了驗(yàn)證這一結(jié)果,圖7給出了不同時(shí)刻平均濃度C=(具體數(shù)值見(jiàn)表2)等值線與靜水壓力拉壓分界線 σH=0 的結(jié)果.可以發(fā)現(xiàn),二者完全重合.

圖7 電極顆粒的彈性分析結(jié)果: 靜水壓力(左)和平均濃度(右)等值線Fig.7 Elastic analysis results of the storage particle: Isolines of the hydrostatic stress (left) and concentration (right)

表2 不同時(shí)刻圓柱體電極顆粒的截面鋰離子平均濃度Table 2 Average concentration of lithium ions in the cylindrical storage particle at various time instants

3.1.2 含初始缺陷的電極顆粒

已有研究表明[13],對(duì)于含初始裂縫的電極顆粒,其破壞模式和破壞程度受幾何尺寸(如半徑r0)、初始裂縫長(zhǎng)度a0、充放電速率c和最大鋰離子濃度Cmax等因素的影響.為驗(yàn)證化學(xué)–力學(xué)耦合相場(chǎng)內(nèi)聚裂縫模型對(duì)于電極顆粒力學(xué)失效行為預(yù)測(cè)的適用性,這里對(duì)含初始缺陷的電極顆粒在脫鋰過(guò)程中的裂縫演化過(guò)程進(jìn)行分析,并考慮三種情況.

(1)初始裂縫不擴(kuò)展.考慮半徑為r0=3 μm 的電極顆粒,其充放電速率c=5,時(shí)間為t=1400 s.對(duì)于初始裂縫長(zhǎng)度為a0=0.1 μm 和a0=1 μm 的兩種電極顆粒,最終時(shí)刻的相場(chǎng)分布和鋰離子濃度分析結(jié)果如圖8 所示.

圖8 含初始缺陷電極顆粒的裂縫不擴(kuò)展模式: 裂縫相場(chǎng)和離子濃度云圖(左 a0=0.1 μm,右 a0=1.0 μm)Fig.8 Pre-notched particles with no cracking growth: Profiles of the crack phase-field and concentration field (left: a0=0.1 μm,right: a0=1.0 μm)

可以看出,對(duì)于有初始裂縫的電極顆粒,雖然應(yīng)力集中導(dǎo)致兩種情況下裂縫尖端的最大主應(yīng)力均超過(guò)材料破壞強(qiáng)度(如圖9 所示),且初始裂縫長(zhǎng)度越大應(yīng)力集中越明顯,但初始裂縫是否擴(kuò)展由能量準(zhǔn)則加以判斷,由于兩種情況下電極顆粒半徑較小且充放電速率較低,直至鋰離子濃度最終近乎為0,流入裂縫尖端的應(yīng)變能仍然不足以驅(qū)動(dòng)初始裂縫進(jìn)一步擴(kuò)展,可以認(rèn)為此時(shí)電極顆粒處于安全狀態(tài)[13].

圖9 含初始缺陷電極顆粒的裂縫不擴(kuò)展模式: 最大主應(yīng)力Fig.9 Pre-notched particles with no cracking growth: Evolution of the maximum principal stress

圖10 給出了電極顆粒內(nèi)部的最大濃度與最小濃度之差 ΔCmax以及平均濃度(通常被稱(chēng)為電極充放電狀態(tài)) SOC 隨時(shí)間的變化曲線.可以看出,初始缺陷長(zhǎng)度對(duì)電極顆粒鋰離子平均濃度的影響不大,但對(duì)最大濃度與最小濃度之差 ΔCmax產(chǎn)生明顯影響.這是因?yàn)? 脫鋰過(guò)程中,濃度梯度 -?C與靜水壓力梯度 ?σH的方向基本一致(如圖11 所示),即靜水壓力對(duì)離子擴(kuò)散起到正向加速的作用,且越靠近邊緣處,擴(kuò)散速率越大.對(duì)于a0=1 μm 的電極顆粒,其應(yīng)力水平更高,電極顆粒表面處離子擴(kuò)散的加速作用更加明顯,由此導(dǎo)致了更大的內(nèi)外濃度差.

圖10 含初始缺陷電極顆粒的裂縫不擴(kuò)展模式Fig.10 Pre-notched particles with no cracking growth

圖11 含初始缺陷電極顆粒的裂縫不擴(kuò)展模式Fig.11 Pre-notched particles with no cracking growth

(2)單裂縫擴(kuò)展.考慮半徑為r0=3 μm、初始裂縫長(zhǎng)度a0=1 μm 的電極顆粒,充放電速率增大至c=10,時(shí)間同樣為t=1400 s.數(shù)值模擬給出的裂縫演化如圖12 所示.

圖12 含初始缺陷電極顆粒的單裂縫擴(kuò)展模式:不同時(shí)刻的裂縫相場(chǎng)云圖Fig.12 Pre-notched storage particles with a single cracking growth:Profiles of the crack phase-field at various time instants

由通量邊界條件 (31)可知: 提高邊界鋰離子的脫出速率會(huì)在電極顆粒表面產(chǎn)生更大的收縮變形和更大的環(huán)向拉應(yīng)力,導(dǎo)致初始裂縫擴(kuò)展;由于電極顆粒內(nèi)部存在壓應(yīng)力區(qū)域,裂縫擴(kuò)展將最終停止(對(duì)于本算例,裂縫擴(kuò)展停止于時(shí)刻t≈177 s,最終裂縫長(zhǎng)度(不包括初始長(zhǎng)度)約為0.55 μm 左右).

(3)多裂縫擴(kuò)展.基于前述分析可知,脫鋰過(guò)程中,隨著裂縫逐漸擴(kuò)展,受內(nèi)部壓應(yīng)力區(qū)影響裂縫尖端處拉應(yīng)力逐漸減小直至裂縫停止發(fā)展;與此同時(shí),由于持續(xù)發(fā)生收縮變形,電極顆粒表面拉應(yīng)力不斷增大.因此,可以預(yù)見(jiàn): 表面拉應(yīng)力增大到一定水平后,可能會(huì)在初始裂縫以外的其他位置產(chǎn)生多條向內(nèi)發(fā)展的附加裂縫.

為了驗(yàn)證該推測(cè),考慮半徑為r0=10 μm、初始裂縫長(zhǎng)度a0=1 μm 的電極顆粒,充放電速率c=10,分析時(shí)間t=50 s.數(shù)值模擬給出的裂縫演化過(guò)程如圖13 所示.

圖13 含初始缺陷電極顆粒的多裂縫擴(kuò)展模式:不同時(shí)刻的裂縫相場(chǎng)云圖Fig.13 Pre-notched particles with multiple cracking growths: Profiles of the crack phase-field at various instants

可以看出,隨著初始裂縫停止擴(kuò)展,電極顆粒表面其他位置處出現(xiàn)新的附加裂縫.其原因,從如圖14所示的裂縫尖端和電極顆粒表面的最大主應(yīng)力隨時(shí)間的變化曲線可以看出:

圖14 含初始缺陷電極顆粒的多裂縫擴(kuò)展模式: 裂縫尖端和電極表面的最大主應(yīng)力隨時(shí)間的變化曲線Fig.14 Pre-notched particles with multiple cracking growths: Evolution of the principal stresses at the crack tip and surface

①時(shí)間段t∈[0 s,9 s] : 初始裂縫尖端和電極顆粒表面的最大主應(yīng)力均呈增長(zhǎng)趨勢(shì),且受應(yīng)力集中影響,裂縫尖端應(yīng)力增長(zhǎng)更快;

② 時(shí)間段t∈[9 s,15 s] : 初始裂縫向內(nèi)部受壓區(qū)發(fā)展,裂縫尖端主拉應(yīng)力水平下降;與此同時(shí),電極表面的拉應(yīng)力繼續(xù)緩慢增加;

③時(shí)間段t∈[15 s,23 s] : 表面鋰離子濃度進(jìn)一步降低,裂縫尖端的最大主應(yīng)力低于材料抗拉強(qiáng)度f(wàn)t=1.3GPa,初始裂縫停止發(fā)展;相應(yīng)地,顆粒表面的主拉應(yīng)力逐漸趨于恒定,但略高于材料強(qiáng)度,電極表面開(kāi)始出現(xiàn)附加裂縫.

3.1.3 無(wú)初始缺陷的電極顆粒

3.1.2 節(jié)分析的電極顆粒均預(yù)設(shè)了一定長(zhǎng)度的初始裂縫,主要考慮脫鋰過(guò)程中初始裂縫是否會(huì)發(fā)生進(jìn)一步擴(kuò)展.本節(jié)進(jìn)一步分析無(wú)初始缺陷的電極顆粒在脫鋰過(guò)程中的力學(xué)失效行為.

(1)安全狀態(tài).考慮半徑為r0=3 μm 的電極顆粒,充放電速率為c=60,分析時(shí)間為t=60 s.數(shù)值模擬給出的裂縫相場(chǎng)演化如圖15 所示.

圖15 無(wú)初始缺陷電極顆粒的安全模式: 裂縫相場(chǎng)演化云圖Fig.15 Intact particles with no cracking grow: Profiles of the crack phase-field at various time instants

可以看出,盡管充放電速率c遠(yuǎn)大于同半徑含初始裂縫的電極顆粒,裂縫相場(chǎng)在t=10 s 左右達(dá)到最大值d≈0.05 后停止增長(zhǎng),此時(shí)電極顆粒表面呈現(xiàn)均勻分布的微裂縫,未發(fā)生局部化并形成向內(nèi)發(fā)展的宏觀裂縫,故可認(rèn)為電極顆粒處于安全狀態(tài).這也說(shuō)明,初始裂縫對(duì)電極顆粒的力學(xué)失效有較大影響.

(2)多裂縫模式.由于電極顆粒表面拉應(yīng)力均勻分布,可以預(yù)期: 隨著電極顆粒尺寸或(和)充放電速率增大,顆粒表面將會(huì)出現(xiàn)多條裂縫.

為了避免對(duì)稱(chēng)模型對(duì)模擬結(jié)果的影響,考慮半徑r0=15 μm 的完整圓形截面模型,充放電速率為c=10,分析時(shí)間為t=200 s.模擬給出的裂縫演化過(guò)程如圖16 所示.

圖16 無(wú)初始缺陷電極顆粒的多裂縫模式: 裂縫相場(chǎng)演化云圖Fig.16 Intact particles with multiple cracking grows: Profiles of the crack phase-field at various time instants

從圖中可以看出: 在鋰離子脫出過(guò)程中,開(kāi)始階段裂縫表面均勻分布的裂縫相場(chǎng)出現(xiàn)局部化,并呈大致對(duì)稱(chēng)形式發(fā)展;擴(kuò)展過(guò)程中,由于內(nèi)部受壓區(qū)的存在,裂縫會(huì)在受壓區(qū)外邊緣沿環(huán)向偏轉(zhuǎn).

為了進(jìn)一步考慮充放電速率的影響,對(duì)同樣半徑的電極顆粒開(kāi)展不同充放電速率下的電極顆粒破壞模式分析.數(shù)值模擬給出的最終裂縫相場(chǎng)分布如圖17所示.

圖17 無(wú)初始缺陷電極顆粒的多裂縫模式: 不同充放電速率下的最終裂縫相場(chǎng)云圖Fig.17 Intact particles with multiple cracking grows: Profiles of the crack phase-field for various charging rates

圖17 無(wú)初始缺陷電極顆粒的多裂縫模式: 不同充放電速率下的最終裂縫相場(chǎng)云圖(續(xù))Fig.17 Intact particles with multiple cracking grows: Profiles of the crack phase-field for various charging rates (continued)

由此可見(jiàn),同等尺寸下,充放電速率越小,出現(xiàn)的總裂縫數(shù)量越多,且每條裂縫幾乎等間距均勻分布于電極顆粒表面,但其最終的裂縫發(fā)展長(zhǎng)度較小;反之則相反.需要指出,當(dāng)充放電速率較高時(shí),每條裂縫為了達(dá)到更大的長(zhǎng)度,裂縫受到內(nèi)部壓應(yīng)力的阻止會(huì)沿環(huán)向偏轉(zhuǎn).

3.1.4 形狀效應(yīng)分析

在彈性階段,已有解析解和數(shù)值結(jié)果表明[29,43],鋰離子擴(kuò)散引起的應(yīng)力不僅與電極顆粒的尺寸大小和充放電速率有關(guān),還受到其截面形狀的影響: 相比截面高寬比較大(如橢圓形)的電極顆粒,高寬比較小(如圓形截面)的電極顆粒會(huì)產(chǎn)生更大的應(yīng)力,由此可能出現(xiàn)更為嚴(yán)重的破壞.

為驗(yàn)證上述結(jié)論,考慮相同面積和邊界鋰離子流量條件下不同高寬比電極顆粒的開(kāi)裂行為.模擬中,電極顆粒面積與半徑r0=5 μm 的圓形電極顆粒相同,式 (31) 中的邊界鋰離子通量均為J0=2.1×10-4mol/(m2·s),脫鋰時(shí)間設(shè)定為100 s.

數(shù)值模擬給出的最終裂縫相場(chǎng)云圖如圖18所示.

圖18 截面面積相同、高寬比不同的電極顆粒: 最終裂縫相場(chǎng)云圖Fig.18 Damage profiles in storage particles with identical section area but of various aspect ratios

相應(yīng)地裂縫長(zhǎng)度隨時(shí)間發(fā)展曲線如圖19 所示.可以看出,同等截面面積和鋰離子邊界條件下,截面高寬比小的電極顆粒更早出現(xiàn)了局部化裂縫,且在相同擴(kuò)散時(shí)間內(nèi)裂縫總長(zhǎng)度更長(zhǎng),因此更容易發(fā)生、更為嚴(yán)重的力學(xué)失效,這與線彈性應(yīng)力分析結(jié)果[29]一致.

圖19 截面面積相同、高寬比不同的電極顆粒: 裂縫長(zhǎng)度演化曲線Fig.19 Crack lengths in storage particles of identical section area but of various aspect ratios

3.2 球體電極顆粒: 三維分析

接下來(lái),考慮球體電極顆粒,其半徑同樣記為r0,相應(yīng)有V/S=r0/3.由于必須采用三維有限元分析,計(jì)算量大大增加,這里僅給出無(wú)初始缺陷電極顆粒的分析結(jié)果.此時(shí),可以采用 1/8 球體進(jìn)行建模,網(wǎng)格劃分采用八節(jié)點(diǎn)實(shí)體單元.

考慮半徑r0=5 μm 的球形電極顆粒,充放電速率為c=20.脫鋰過(guò)程中電極顆粒內(nèi)部的鋰離子濃度分布如圖20 所示.與二維情況一致,鋰離子濃度分布由內(nèi)向外單調(diào)遞減.

圖20 球體電極顆粒三維分析: 鋰離子濃度分布示意圖Fig.20 Spherical storage particles: Concentration contour of Li-ion

3.2.1 彈性分析結(jié)果

在彈性階段,電極顆粒的主應(yīng)力分布及方向如圖21 所示.

圖21 球體電極顆粒三維分析: 彈性階段的主應(yīng)力云圖Fig.21 Spherical storage particles: Contours of principal stresses in the elastic stage

可以看出,球體電極顆粒三維彈性分析結(jié)果與圓柱體二維分析類(lèi)似,應(yīng)力狀態(tài)均為外部受拉、內(nèi)部受壓,第一、二主應(yīng)力分別沿經(jīng)線和緯線方向,其大小相同且大于徑向的第三主應(yīng)力.這一結(jié)果與附錄給出的解析解一致.

3.2.2 損傷破壞全過(guò)程分析

接下來(lái),利用相場(chǎng)內(nèi)聚裂縫模型分析球體電極顆粒的損傷破壞全過(guò)程.由于前期大量二維和三維分析結(jié)果表明[33,44-45],相場(chǎng)內(nèi)聚裂縫模型中尺度參數(shù)b的取值僅影響裂縫帶的寬度,對(duì)破壞模式和整體計(jì)算結(jié)果幾乎沒(méi)有影響.因此,為減少計(jì)算量,三維分析中取用了稍大的相場(chǎng)尺度參數(shù)b=0.14 μm.

三維數(shù)值模擬給出的電極顆粒最終裂縫相場(chǎng)分布如圖22 所示.

圖22 球體電極顆粒三維分析: 裂縫相場(chǎng)分布云圖Fig.22 Spherical storage particles: Profiles of the crack phase-field

結(jié)合圖22 中給出的主應(yīng)力分布,最終破壞時(shí)球體電極顆粒表面的裂縫帶大致可以歸為兩類(lèi): (1) 靠近上、下兩極處沿緯度方向的兩條環(huán)向裂縫,受第一主拉應(yīng)力驅(qū)動(dòng);(2)沿經(jīng)線方向的豎直裂縫,受第二主拉應(yīng)力驅(qū)動(dòng).由于第一、二主拉應(yīng)力大小相等而方向正交,故兩類(lèi)裂縫幾乎同時(shí)出現(xiàn),最終分布也近似垂直.

為進(jìn)一步研究球體電極顆粒內(nèi)部的裂縫演化過(guò)程,沿X-Z平面切割球體,其剖面上的裂縫相場(chǎng)演化過(guò)程如圖23 所示.有趣的是,該剖面上的裂縫相場(chǎng)演化與二維圓柱體電極顆粒(無(wú)初始缺陷的完好電極)基本相同.究其原因,對(duì)于X-Z剖平面,與之垂直的第二主拉應(yīng)力對(duì)該平面上的裂縫演化不起貢獻(xiàn),即裂縫發(fā)展僅由第一主拉應(yīng)力驅(qū)動(dòng),因此其裂縫演化情況與二維分析結(jié)果一致.

圖23 球體電極顆粒三維分析: X-Z 平面裂縫相場(chǎng)演化Fig.23 Spherical storage particles: Evolution of the crack phase-field on the X-Z plane

4 結(jié)論與展望

在統(tǒng)一相場(chǎng)理論框架內(nèi),本文建立了化學(xué)擴(kuò)散-力學(xué)變形耦合的相場(chǎng)內(nèi)聚裂縫模型,發(fā)展了相應(yīng)的多場(chǎng)有限元數(shù)值實(shí)現(xiàn)算法,并應(yīng)用于鋰電池電極顆粒裂縫演化導(dǎo)致的力學(xué)失效分析.

與純力學(xué)作用下的相場(chǎng)內(nèi)聚裂縫模型一樣,上述模型同時(shí)融入了基于強(qiáng)度的裂縫起裂準(zhǔn)則、基于能量的裂縫擴(kuò)展準(zhǔn)則以及基于變分原理的裂縫擴(kuò)展路徑判據(jù),因此能夠描述鋰離子脫出過(guò)程中電極顆粒的復(fù)雜裂縫演化過(guò)程.

數(shù)值算例表明,所提出的化學(xué)-力學(xué)耦合相場(chǎng)內(nèi)聚裂縫模型不僅能夠描述帶初始缺陷電極顆粒的開(kāi)裂行為,而且同樣適用于無(wú)初始缺陷的完好電極顆粒.由于數(shù)值模擬結(jié)果不受裂縫尺度和網(wǎng)格大小的影響,該模型能夠應(yīng)用于鋰離子擴(kuò)散過(guò)程中電極顆粒損傷破壞全過(guò)程的二維和三維分析,對(duì)商業(yè)鋰電池電極的優(yōu)化設(shè)計(jì)具有一定的指導(dǎo)意義.

需要指出的是,本文忽略了鋰離子濃度變化對(duì)電極顆粒材料力學(xué)參數(shù)(如彈性模量等)的影響,且尚未考慮溫度場(chǎng)對(duì)濃度分布和裂縫演化的耦合作用.這些問(wèn)題將在后續(xù)研究中進(jìn)一步加以考慮.

附錄

為完整起見(jiàn),這里給出圓柱體和球體兩種形狀的電極顆粒在彈性階段的應(yīng)力–離子濃度關(guān)系解析解.

(1) 圓柱體電極顆粒

采用極坐標(biāo) (r,φ) 表示圓柱體電極顆粒的平面應(yīng)力狀態(tài).對(duì)于截面形狀為圓形的電極顆粒,任意一點(diǎn)的剪應(yīng)力 τrφ=0,應(yīng)力狀態(tài)由徑向應(yīng)力 σr和環(huán)向應(yīng)力 σφ確定,且其應(yīng)力狀態(tài)不隨角度 φ 改變,僅為徑向距離r的函數(shù).相應(yīng)地,變形協(xié)調(diào)方程和線彈性應(yīng)力–應(yīng)變關(guān)系為

或表示為如下等效形式

其中,E0和 ν0分別為材料的彈性模量和泊松比.將式 (A2)代入靜力平衡方程

式中的積分常數(shù)K1和K2由以下邊界條件確定:

①電極顆粒表面不受任何面力作用,徑向應(yīng)力為0,即σr|(r=r0)=0;

② 電極顆粒中心點(diǎn)r=0 處應(yīng)力取極值(壓應(yīng)力最大),即

相應(yīng)地,K1和K2確定為

于是,徑向應(yīng)力 σr、環(huán)向應(yīng)力 σφ 和靜水壓力 σH分別表示為

即靜水應(yīng)力 σH的拉壓分界線 σH(r)=0 與電極顆粒的整體平均濃度等值線C(r)=Cˉ(r0) 重合.同時(shí),對(duì)于脫鋰過(guò)程,濃度分布從中心到外邊界單調(diào)遞減,故有 σφ≥σr,即最大主應(yīng)力為環(huán)向應(yīng)力.

(2) 球體電極顆粒

采用球坐標(biāo) (r,θ,φ) 表示球體電極顆粒的應(yīng)力狀態(tài).此時(shí),根據(jù)類(lèi)似圓柱形顆粒的應(yīng)力分析方法,可得出其極坐標(biāo)下的應(yīng)力解析解.由于幾何中心對(duì)稱(chēng)性,剪應(yīng)力 τrθ和 τrφ為零,應(yīng)力狀態(tài)由一個(gè)徑向應(yīng)力 σr和兩個(gè)環(huán)向應(yīng)力 σθ和 σφ完全確定且兩環(huán)向應(yīng)力在數(shù)值上相等.此外,對(duì)于球體幾何形狀,任意一點(diǎn)的應(yīng)力狀態(tài)不隨轉(zhuǎn)角 θ 和 φ 發(fā)生改變,僅為徑距r的函數(shù).

球體電極顆粒應(yīng)力解析解推導(dǎo)過(guò)程與二維平面應(yīng)力問(wèn)題類(lèi)似,最終其應(yīng)力–離子濃度關(guān)系表示為

即靜水應(yīng)力 σH的拉壓分界線 σH(r)=0 同樣與電極顆粒的整體平均濃度等值線重合,且最大主應(yīng)力同樣為環(huán)向應(yīng)力即 σφ≥σr.