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

謝奕展 程夕明

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

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

1 引言

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

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

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

為了以盡量少的時間準確求解ESPM 液相擴散過程，假設電化學反應僅發生在電解質與集電極界面處，從而僅需求解界面處鋰離子濃度.為了高效、準確地獲得恒流工況和動態工況下的界面鋰離子濃度，一種基于液相擴散方程解析解的有限離散卷積方法被提出.

2 模型建立

2.1 問題描述

在改進單粒子模型中，考慮液相的影響，并假設正負極及隔膜的電導率、擴散系數等電化學性質相同，且保持不變，那么液相的質量守恒表達為

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

電勢守恒方程為

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

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

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

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

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

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

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

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

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

由(1)式和(3)式聯立，得

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

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

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

2.2 描述方程及其解析解

2.2.1 恒流工況

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

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

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

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

則液相兩側濃度的拉普拉斯變換為

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

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

并令

由于

且{fN(0)}是一單調遞增數列，所以

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

因此，(17)式可重寫為

令

那么，(22)式可表達為

2.2.2 動態工況

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

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

由(24)式得系統的單位脈沖響應為

則液相L側的濃度為

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

記

則(28)式可重寫為

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

2.3 動態工況求解

2.3.1 離散化

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

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

而

記另一收斂和函數:

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

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

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

2.3.2 截 斷

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

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

表1 模型參數*Table 1.Model Parameters*.

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

3 仿真與結果

3.1 實驗仿真說明

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

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

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

3.2 誤差分析

3.2.1 恒流工況

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

圖3 恒流工況仿真結果 (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，且均發生在初始時刻.隨著時間增大，誤差逐漸變小為零.注意到最大誤差之比近似為對應放電倍率之比，其原因可能在于Comsol 軟件處理階躍輸入時使用了一定的平滑過渡處理，使得Comsol 的計算結果偏大，而這一平滑處理導致誤差與倍率成正比.但是，隨著時間的推移，輸入逐漸變為恒定值，使得Comsol 的計算結果與理論解結果相同.與專業軟件Comsol 的有限元方法相比，在恒流工況下，所提解法是一種解析解，僅需(24)式中g(t)函數以及相關系統參數，計算簡便.

3.2.2 動態工況

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

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

對比圖4(b)、圖5(b)和圖6(b)，當N0=200時，各動態工況的的本算法鋰離子濃度曲線AnaN0(200)DST，AnaN0(200)NEDC，AnaN0(200)WLTC與其對應Comsol 仿真的濃度曲線eT(1×10—5)DST，eT(1×10—5)NEDC，eT(1×10—5)WLTC 差別較大，而當N0=∞時，即h(t)無截斷時，本算法解與對應的Comsol 仿真結果一致.同樣，由圖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 時截斷將導致誤差較大.分析圖4(b)、圖5(b)和 圖6(b)，與Comsol 仿真結果eT(1×10—5)DST，eT(1×10—5)NEDC 和eT(1×10—5)WLTC相比，所提方法在N0=200 時的濃度曲線與Comsol仿真結果eT(1×10—5)DST，eT(1×10—5)NEDC，eT(1×10—5)WLTC 相差最大，N0=600 時次之，而N0=1000 和N0=∞時均與Comsol 仿真結果高度一致.出現上述誤差現象是由圖1 的h(t)在不同截斷時長截斷所產生.顯然，N0越小，截斷時長T越小，h(T)越小，導致(38)式截斷為(40)式所忽略的項:

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

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

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

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

表2 動態工況誤差Table 2.Errors under dynamic profiles.

3.3 計算效率分析

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

表3 列出了5 組實驗的誤差RRMSE(%)值和計算效率指標.其中，隨著N0的變大，卷積算法的求解時間有所增大，這與(40)式的計算規律相符合.在(40)式中，當N0增大時，所需計算的kN+1-n(hn-hn-1)項增多，導致求解時間有所增長，但相對求解時間都不超過1/1000.在N0=600 時，卷積算法的結果AnaN0(600)WLTC 與Comsol 結果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 結論

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

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

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