一種求解鋰離子電池單粒子模型液相擴(kuò)散方程的新方法*

謝奕展 程夕明

(北京理工大學(xué)機械與車輛學(xué)院，電動車輛國家工程實驗室，北京 100081)

電解液中的鋰離子濃度表達(dá)是鋰離子電池電化學(xué)模型求解的基本任務(wù)之一.為了平衡單粒子模型的液相動態(tài)性能和計算效率，假設(shè)反應(yīng)僅發(fā)生在集電極和電解質(zhì)界面上，為此，提出一種基于液相擴(kuò)散方程無窮級數(shù)解析解的界面濃度求解新方法.在恒流工況下，利用數(shù)列單調(diào)收斂準(zhǔn)則將解析解轉(zhuǎn)化為一個收斂和函數(shù).在動態(tài)工況下，將該解析解簡化為輸入與和函數(shù)的無限離散卷積.利用和函數(shù)隨時間單調(diào)衰減并收斂至零的特性對其進(jìn)行截斷，從而得到有限離散卷積求解算法.對比專業(yè)有限元分析軟件，該方法在恒流工況和動態(tài)工況下均能以較少的計算時間獲得相當(dāng)好的精度.而且，該方法僅有一個配置參數(shù).因此，所提方法將有效減小應(yīng)用于實時電池管理系統(tǒng)上的電化學(xué)模型計算負(fù)擔(dān).

1 引言

鋰離子電池具有能量密度高、功率密度高、循環(huán)壽命長等特點[1]，目前是純電動汽車的主要動力能源裝置.整車的高效安全運行依賴鋰離子電池狀態(tài)，如電壓、電流、內(nèi)部溫度、荷電狀態(tài)、健康狀態(tài)等的有效監(jiān)測和管理[2].因此，十分必要建立準(zhǔn)確高效的鋰離子電池模型以可靠估計電池狀態(tài).

目前，鋰離子電池建模方法主要有機器學(xué)習(xí)方法、等效電路模型和電化學(xué)模型三大類.機器學(xué)習(xí)方法通過支持向量機、相關(guān)向量機、高斯過程、神經(jīng)網(wǎng)絡(luò)等方法自適應(yīng)習(xí)得特征參數(shù)與狀態(tài)量的函數(shù)關(guān)系[3-6].機器學(xué)習(xí)方法不需要了解電池內(nèi)部復(fù)雜的電化學(xué)過程，具有實施簡單、門檻低的優(yōu)點.但是，機器學(xué)習(xí)方法嚴(yán)重依賴樣本數(shù)據(jù)，其模型參數(shù)缺少明確的物理意義，普適性差.等效電路模型將電池視作多個電阻電容串并聯(lián)的系統(tǒng)，并利用實驗數(shù)據(jù)辨識模型參數(shù)[7-9].等效電路模型具有參數(shù)少、模型簡單、易于實時計算的優(yōu)點，但是其參數(shù)同樣缺少具體的物理意義，且難以寬范圍適應(yīng)電池的時變非線性行為.電化學(xué)模型通過數(shù)理方程描述電池內(nèi)部工作機理，具有明確的物理意義.為此，一些基于偽二維模型(pesudo-two-dimensional model，P2D)的電化學(xué)降階模型被應(yīng)用于鋰離子電池建模中[10，11].其中，由于計算量較小，單粒子模型(single particle model，SPM)被廣泛應(yīng)用[12，13].然而，SPM由于忽略液相擴(kuò)散過程而導(dǎo)致其在高放電倍率下誤差較大.為此，許多改進(jìn)單粒子模型(enhanced single particle model，ESPM)考慮了液相擴(kuò)散過程以提高模型精度[14，15].

ESPM 實時應(yīng)用的難點在于固相擴(kuò)散方程和液相物質(zhì)傳輸方程的快速而高精度的計算.其中，基于球形粒子的電極固相擴(kuò)散方程的實時求解問題已經(jīng)解決，至少可采用多項式近似、Padé近似和基于解析解的離散卷積等三種方法高效而可靠求解[16].針對ESPM 的液相擴(kuò)散方程的求解問題，主要的求解方法有解析法[17]、微元法、傳遞函數(shù)近似法[18，19]和多項式近似法[14，15，20].盡管解析法能夠得到準(zhǔn)確的結(jié)果，但是該方法的動態(tài)工況表達(dá)式十分復(fù)雜[17].微元法主要有有限元方法、有限體積法和有限差分法，并已被廣泛應(yīng)用于商業(yè)軟件中，然而它們的參數(shù)眾多且計算量大.為了降低計算負(fù)擔(dān)，整數(shù)階傳遞函數(shù)近似和多項式近似被廣泛使用，然而它們在一般工況下的精度較低，并且傳遞函數(shù)和多項式形式的確定方法缺少一般性[21].

為了以盡量少的時間準(zhǔn)確求解ESPM 液相擴(kuò)散過程，假設(shè)電化學(xué)反應(yīng)僅發(fā)生在電解質(zhì)與集電極界面處，從而僅需求解界面處鋰離子濃度.為了高效、準(zhǔn)確地獲得恒流工況和動態(tài)工況下的界面鋰離子濃度，一種基于液相擴(kuò)散方程解析解的有限離散卷積方法被提出.

2 模型建立

2.1 問題描述

在改進(jìn)單粒子模型中，考慮液相的影響，并假設(shè)正負(fù)極及隔膜的電導(dǎo)率、擴(kuò)散系數(shù)等電化學(xué)性質(zhì)相同，且保持不變，那么液相的質(zhì)量守恒表達(dá)為

式中，j(t，x)為鋰離子在液相厚度方向x處在t時刻的通量密度(mol·m—2·s—1);D為擴(kuò)散系數(shù)(m2/s);c(x，t)為液相厚度方向x處在t時刻的濃度(mol/m3);i(t，x)為液相厚度方向x處在t時刻的電流密度大小(A/m2);t+為離子遷移數(shù)，量綱為1;F是法拉第常數(shù);L為液相厚度R的一半(m).

電勢守恒方程為

式中，σ為液相電導(dǎo)率(S/m);φ(x，t)為液相厚度方向x處在t時刻的電勢(V);Rg為氣體常數(shù);T為溫度(K);在改進(jìn)單粒子模型中，如圖1 液相模型所示，假設(shè)固相與液相鋰離子僅在液相邊界處，即x=—L，L處進(jìn)入，則邊界條件為

圖1 液相模型Fig.1.The electrolyte model.

式中，j(t)為正負(fù)極處固相與液相的邊界鋰離子交換通量密度(mol·m—2·s—1).

假設(shè)初始時刻濃度處處分布均勻，故初值條件為

規(guī)定—L處電勢為零點，故另一邊界條件為

由于考慮液相的單粒子模型中固相與液相僅在邊界處發(fā)生反應(yīng)，根據(jù)Bulter-Volmer 方程，液相僅需要求取液相兩側(cè)的鋰離子濃度和電勢，也就是c(t，L)，c(t，—L)，φ(t，L)和φ(t，—L).

聯(lián)立(1)式、(2)式、(3)式、(5)式和(6)式，得

由(8)式可求解得到c(t，L)和c(t，—L).

將(3)式代入(4)式，并對(4)式進(jìn)行積分，得

由(1)式和(3)式聯(lián)立，得

將(7)式和(10)式代入(9)式中，有

由(11)式可知，φ(t，L)可由c(t，L)和c(t，—L)直接求得.

綜上，改進(jìn)單粒子模型液相部分僅需求解式(8)，從而得到液相兩側(cè)的濃度c(t，L)和c(t，—L).

2.2 描述方程及其解析解

2.2.1 恒流工況

在(8)式中，遷移系數(shù)t+能被整合進(jìn)擴(kuò)散系數(shù)D，等效為一個新的擴(kuò)散系數(shù)，因此將t+視為0 并不影響求解(8)式方法的一般性.假設(shè)初始濃度為0，那么在恒流工況下，鋰離子的液相擴(kuò)散過程可由偏微分方程及其初始條件、邊界條件描述:

式中，jconst為液相兩側(cè)界面通量密度(mol·m—2·s—1).

對(12)式進(jìn)行拉普拉斯變換，得

式中，C(x)是c(t，x)的拉普拉斯變換.求解(13)式，得

則液相兩側(cè)濃度的拉普拉斯變換為

顯然，C(L)和C(—L)互為相反數(shù)，因此在下文內(nèi)容中將僅考慮C(L).對(16)式進(jìn)行拉普拉斯逆變換，得L側(cè)的濃度為

由于(17)式是一無窮級數(shù)，需要進(jìn)行簡化以方便求解.令收斂和函數(shù):

并令

由于

且{fN(0)}是一單調(diào)遞增數(shù)列，所以

又因為fN(t) ＞ 0 且隨t單調(diào)衰減，所以對任意t＞ 0，數(shù)列{fN(t)}均單調(diào)遞增且有上下界.從而，{fN(t)}存在極限，即(18)式的f(t)存在且是一有界函數(shù).

因此，(17)式可重寫為

令

那么，(22)式可表達(dá)為

2.2.2 動態(tài)工況

在動態(tài)工況下，(12)式改寫為

式中j(t)為液相兩側(cè)端點時變通量密度(mol·m—2·s—1).

由(24)式得系統(tǒng)的單位脈沖響應(yīng)為

則液相L側(cè)的濃度為

對(27)式進(jìn)行分部積分，得

記

則(28)式可重寫為

顯然，(30)式的求解主要在于(29)式中y(t)的求解.

2.3 動態(tài)工況求解

2.3.1 離散化

(29)式的y(t)是動態(tài)工況濃度求解(30)式的難點.為此，將區(qū)間[0，t]等分為N個區(qū)間，第n個區(qū)間為[tn-1，tn]，tn=ndt，每個區(qū)間大小為dt=t/N.目前，世界輕型汽車測試規(guī)程(worldwide harmonized light vehicles test cycle，WLTC)、新標(biāo)歐洲循環(huán)測試(new European driving cycle，NEDC)工況等均以1 s 為一個采樣點，因此選取dt=1 s，并假設(shè)節(jié)點之間的值由線性插值所得，并簡記j(tn)=jn，則有

式中，kn是j(t)在區(qū)間[tn-1，tn]上的斜率.從而，(29)式可改寫為

而

記另一收斂和函數(shù):

與f(t)的分析類似，h(t)極限存在且單調(diào)有界.簡記hn=h(tn)，有

將(36)式代入(33)式中，有

將(37)式代入(30)式中，并簡記g(tn)=gn，整理得動態(tài)工況下濃度的離散卷積形式解析解為

2.3.2 截 斷

顯然(38)式的計算量隨著t的增大而線性增長，使得其在t很大時不具備可行性，需要對其進(jìn)行簡化.注意到隨h(t)隨t單調(diào)衰減，且有

根據(jù)表1 參數(shù)繪制得h(t)的函數(shù)圖像如圖2所示.顯然，當(dāng)t大于某個值，稱為截斷時長T，且假設(shè)此時N0=T/dt，可認(rèn)為n≥N0時，hn—hn—1≈0，并考慮初始濃度，由(38)式簡化得濃度的有限離散卷積算法為:

表1 模型參數(shù)*Table 1.Model Parameters*.

圖2 h(t)函數(shù)曲線Fig.2.The curve of h(t).

3 仿真與結(jié)果

3.1 實驗仿真說明

由于放電和充電工況的液相兩側(cè)表面通量的方向相反，因此以放電和充放電混合工況分析求解方法的性能.一般地，純電動汽車動力電池系統(tǒng)的最大放電倍率不超過3 C，為此，采用6 種工況，包括0.25 C，1 C 和3 C 三種恒流放電工況和以最高倍率為3 C 的NEDC，WLTC，DST 三種動態(tài)工況，驗證模型求解方法，并與有限元方法(Comsol軟件)的計算結(jié)果對比，分析兩者的模型誤差和計算時間，其中表1 列示了模型參數(shù).電腦配置為AMD Ryzen 5 4500 U，Radeon Graphics 2.38 GHz 處理器和16.0 GB (15.2 GB 可用)內(nèi)存.

為權(quán)衡仿真結(jié)果收斂情況和計算時間，將Comsol 軟件的容許誤差eT(torelence error)設(shè)置為1×10—5.若無特殊說明，Comsol 軟件的實驗結(jié)果均在eT=1×10—5下取得，并將其視為參考解.本文采用以下2 種無量綱誤差:1) 相對均方根誤差(reference root mean squard error，RRMSE)，為均方根誤差與表1 中的初始濃度之比;2) 最大相對誤差(reference maximum absolute error，RMAE)，為最大絕對誤差與初始濃度之比.

式中，c0為表1中的初始濃度;m為樣本數(shù)目;ui為所提方法計算值;vi為參考解.

3.2 誤差分析

3.2.1 恒流工況

在恒流工況下，電極兩側(cè)濃度表達(dá)式為(24)式.圖3 分別顯示了三個恒流工況0.25 C，1 C 和3 C的實驗結(jié)果.圖3(a)、圖3(c)和圖3(e)顯示了三種恒流倍率的理論解濃度曲線為AnaC(0.25)，AnaC(1)，AnaC(3)及其對應(yīng)的Comosl 仿真所得濃度曲線eT(1×10—5)C(0.25)，eT(1×10—5)C(1)和eT(1×10—5)C(3)，它們高度吻合.

圖3 恒流工況仿真結(jié)果 (a) 0.25 C 濃度;(b) 0.25 C 誤差;(c) 1.0 C 濃度;(d) 1.0 C 誤差;(e) 3.0 C 濃度;(f) 3.0 C 誤差Fig.3.Simulation results of the galvanostatic profiles:(a) 0.25 C concentration;(b) 0.25 C error;(c) 1.0 C concentration;(d) 1.0 C error;(e) 3.0 C concentration;(f) 3.0 C error.

由圖3(b)、圖3(d)和圖3(f)可知，三種倍率的最大誤差依次為—5.5 mol/m3，—22.5 mol/m3和—67.8 mol/m3，且均發(fā)生在初始時刻.隨著時間增大，誤差逐漸變小為零.注意到最大誤差之比近似為對應(yīng)放電倍率之比，其原因可能在于Comsol 軟件處理階躍輸入時使用了一定的平滑過渡處理，使得Comsol 的計算結(jié)果偏大，而這一平滑處理導(dǎo)致誤差與倍率成正比.但是，隨著時間的推移，輸入逐漸變?yōu)楹愣ㄖ?，使得Comsol 的計算結(jié)果與理論解結(jié)果相同.與專業(yè)軟件Comsol 的有限元方法相比，在恒流工況下，所提解法是一種解析解，僅需(24)式中g(shù)(t)函數(shù)以及相關(guān)系統(tǒng)參數(shù)，計算簡便.

3.2.2 動態(tài)工況

在動態(tài)工況下，電極兩側(cè)鋰離子濃度由(40)式計算得到.本節(jié)將在DST，NEDC 和WLTC 三種工況，并分別將(40)式的N0值選取為200，600，1000 和∞(表示h(t)無截斷)，對該算法進(jìn)行精度和求解時間的驗證.由于一個DST 工況持續(xù)時長僅360 s，僅存在N0=200 和N0=∞兩種情況.而NEDC 和WLTC 工況各有四種情況，一共有10 組實驗.圖4、圖5 和圖6 依次顯示了DST，NEDC，WLTC 工況以及相對應(yīng)的濃度曲線和誤差曲線.

圖4 DST 工況仿真結(jié)果 (a) 放電倍率;(b) 濃度;(c) 誤差Fig.4.Simulation results for the DST profile:(a) Profile;(b) concentration;(c) error.

對比圖4(b)、圖5(b)和圖6(b)，當(dāng)N0=200時，各動態(tài)工況的的本算法鋰離子濃度曲線AnaN0(200)DST，AnaN0(200)NEDC，AnaN0(200)WLTC與其對應(yīng)Comsol 仿真的濃度曲線eT(1×10—5)DST，eT(1×10—5)NEDC，eT(1×10—5)WLTC 差別較大，而當(dāng)N0=∞時，即h(t)無截斷時，本算法解與對應(yīng)的Comsol 仿真結(jié)果一致.同樣，由圖4(c)、圖5(c)和圖6(c)可知，所提算法在N0=200 時的誤差曲線Err(Ana-eT)N0(200)DST，Err(Ana-eT)N0(200)NEDC 和 Err(Ana-eT)N0(200)WLTC 比在N0=∞時的誤差曲線波動要劇烈得多，說明h(t)在T=200 s 時截斷將導(dǎo)致誤差較大.分析圖4(b)、圖5(b)和 圖6(b)，與Comsol 仿真結(jié)果eT(1×10—5)DST，eT(1×10—5)NEDC 和eT(1×10—5)WLTC相比，所提方法在N0=200 時的濃度曲線與Comsol仿真結(jié)果eT(1×10—5)DST，eT(1×10—5)NEDC，eT(1×10—5)WLTC 相差最大，N0=600 時次之，而N0=1000 和N0=∞時均與Comsol 仿真結(jié)果高度一致.出現(xiàn)上述誤差現(xiàn)象是由圖1 的h(t)在不同截斷時長截斷所產(chǎn)生.顯然，N0越小，截斷時長T越小，h(T)越小，導(dǎo)致(38)式截斷為(40)式所忽略的項:

圖5 NEDC 工況仿真結(jié)果 (a) 放電倍率;(b) 濃度;(c) 誤差Fig.5.Simulation results for the NEDC profile:(a) Profile;(b) concentration;(c) error.

圖6 WLTC 工況仿真結(jié)果 (a) 放電倍率;(b) 濃度;(c) 誤差Fig.6.Simulation results for the WLTC profile:(a) Profile，(b) concentration，(c) error.

越大，使截斷誤差增大.而當(dāng)N0=1000 和N0=∞時，卷積算法解均與Comsol 仿真結(jié)果高度一致，原因在于此時的h(T)十分接近于零，忽略項(43)式很小，從而解析解誤差與沒截斷時差別不大且十分小.

表2 列出了各動態(tài)工況的誤差統(tǒng)計結(jié)果.其中，在三種動態(tài)工況下，當(dāng)N0=∞，即h(t)無截斷時，卷積算法在DST，NEDC 和WLTC 工況下的相對均方根誤差RRMSE 依次為0.09%，0.04%和0.05%，RMAE 依次為0.61%，0.18%和0.16%.相比于其他工況，DST 工況的求解誤差最大，其原因可能是DST 工況的電流階躍幅度正負(fù)變化最大，從而由Comsol 平滑處理大階躍時導(dǎo)致的誤差增大.當(dāng)h(t)的截斷時長T越小時，誤差越大.在NEDC 和WLTC 工況下，當(dāng)N0=200 時，各動態(tài)工況的濃度誤差RRMSE 是N0=600 的19 倍，N0=600 僅為N0=1000 的3 倍，而N0=1000與N0=∞二者誤差可忽略不計.當(dāng)N0取值較小時，提高N0能有效提高卷積算法精度;當(dāng)N0較大時，提高N0并不能提高精度，說明截斷h(t)能夠在不犧牲求解精度前提下減小算法計算量.上述規(guī)律符合h(t)衰減指數(shù)化下降并逼近于0 的性質(zhì).誤差最大時發(fā)生在WLTC 工況下N0取200 時，為3%，最小時發(fā)生在NEDC 工況下N0=1000 時，僅為0.04%.

表2 動態(tài)工況誤差Table 2.Errors under dynamic profiles.

3.3 計算效率分析

以WLTC 工況為例，比較Comsol 容許誤差設(shè)置為eT=1×10—3，以及所提算法在N0=200，600，1000 和∞共5 組實驗，并將Comsol 容許誤差設(shè)置為eT=1×10—5的濃度作為參考值.同時，將eT=1×10—5的求解時間作為參考值，采取相對求解時間RT(reference time，表示實驗求解時間與eT=1×10—5求解時間之比) 作為衡量求解時間的指標(biāo).當(dāng)容許誤差設(shè)置為eT(1×10—5)時，WLTC 工況的Comsol 求解時間為5.18 s.

表3 列出了5 組實驗的誤差RRMSE(%)值和計算效率指標(biāo).其中，隨著N0的變大，卷積算法的求解時間有所增大，這與(40)式的計算規(guī)律相符合.在(40)式中，當(dāng)N0增大時，所需計算的kN+1-n(hn-hn-1)項增多，導(dǎo)致求解時間有所增長，但相對求解時間都不超過1/1000.在N0=600 時，卷積算法的結(jié)果AnaN0(600)WLTC 與Comsol 結(jié)果eT(1×10—3)WLTC 相比，二者求解精度接近，但所提卷積算法的求解時間僅為Comsol 的1/500.

表3 WLTC 工況下卷積算法與Comsol 比較Table 3.The comparison between the convolution algorithm and COMSOL under the WLTC.

4 結(jié)論

本文提出了一種新的求解鋰離子電池電化學(xué)模型液相擴(kuò)散過程液相兩側(cè)表面濃度的離散卷積算法.該算法求解簡單，參數(shù)較少，僅需配置1 個參數(shù)，模型的求解誤差隨著截斷時長T的增加而減小.

該離散卷積算法是基于擴(kuò)散方程的解析解而提出，具有收斂快速而精度高的優(yōu)點.在DST，NEDC和WLTC 工況下，相比于容許誤差eT=1×10—5的Comsol 模型求解結(jié)果，所提算法取截斷誤差最大時的計算誤差也僅為3%，而取截斷誤差最小時的計算誤差僅為0.04%左右，具有相當(dāng)高的計算精度.在WLTC 工況激勵下，該卷積算法能夠以有限元方法1/500 的求解時間取得與后者相近的精度.

在數(shù)據(jù)存儲成本方面，該離散卷積算法僅需存儲工況數(shù)據(jù)、收斂和函數(shù)g(t)和h(t)的各N0個樣本值，數(shù)據(jù)存儲量少.因此，所提方法將有助于電化學(xué)模型在電池管理系統(tǒng)上的實時應(yīng)用.對于不同的電池液相參數(shù)，特別是時變參數(shù)，算法的適用性需進(jìn)一步研究.在下一階段工作中，將聚焦于該離散卷積算法在鋰離子電池電化學(xué)降階模型上的集成.