吳建營 洪屹峰

* (華南理工大學(xué)亞熱帶建筑科學(xué)國家重點實驗室,廣州 510641)

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

引言

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

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

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

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

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

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

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

雖然統(tǒng)一相場理論也適用于有限變形[28],這里僅考慮小應(yīng)變情況,故采用位移場u(x) 和線性應(yīng)變場 ε(x)=?su(x) 描述固體的變形狀態(tài),這里 ?s(·) 表示對稱梯度算子.伴隨鋰電池充放電過程中電極顆粒表面鋰離子的脫出或嵌入,離子濃度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為二階單位張量.

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

1.1 擴散–變形耦合過程的控制方程

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

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

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

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

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

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

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

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

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

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

1.2 裂縫相場演化的能量原理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖1 鋰電池電極顆粒的擴散-開裂-變形耦合過程Fig.1 Chemo-mechanical coupling in storage particles of LIBs

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

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

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

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

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

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

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

2.1 控制方程的弱形式

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

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

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

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

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

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

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

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

2.3 時間域離散

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

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

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

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

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

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

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

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

2.4 交錯迭代算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖5 脫鋰過程電極顆粒的應(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,分析時間為t=5000 s.圖6 給出了最大主應(yīng)力隨時間的變化曲線.

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

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

在線彈性階段,附錄給出的解析解表明,圓柱體電極顆粒靜水應(yīng)力的拉壓分界線 σH(r)=0 與截面平均濃度等值線C(r)=重合.為了驗證這一結(jié)果,圖7給出了不同時刻平均濃度C=(具體數(shù)值見表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 不同時刻圓柱體電極顆粒的截面鋰離子平均濃度Table 2 Average concentration of lithium ions in the cylindrical storage particle at various time instants

3.1.2 含初始缺陷的電極顆粒

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

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

圖8 含初始缺陷電極顆粒的裂縫不擴展模式: 裂縫相場和離子濃度云圖(左 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)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1.3 無初始缺陷的電極顆粒

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

考慮半徑r0=5 μm 的球形電極顆粒,充放電速率為c=20.脫鋰過程中電極顆粒內(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é)果與圓柱體二維分析類似,應(yīng)力狀態(tài)均為外部受拉、內(nèi)部受壓,第一、二主應(yīng)力分別沿經(jīng)線和緯線方向,其大小相同且大于徑向的第三主應(yīng)力.這一結(jié)果與附錄給出的解析解一致.

3.2.2 損傷破壞全過程分析

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

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

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

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

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

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

4 結(jié)論與展望

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

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

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

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

附錄

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

(1) 圓柱體電極顆粒

采用極坐標 (r,φ) 表示圓柱體電極顆粒的平面應(yīng)力狀態(tài).對于截面形狀為圓形的電極顆粒,任意一點的剪應(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;

② 電極顆粒中心點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) 重合.同時,對于脫鋰過程,濃度分布從中心到外邊界單調(diào)遞減,故有 σφ≥σr,即最大主應(yīng)力為環(huán)向應(yīng)力.

(2) 球體電極顆粒

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

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

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