搜索

x

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

一种求解锂离子电池单粒子模型液相扩散方程的新方法

谢奕展 程夕明

引用本文:
Citation:

一种求解锂离子电池单粒子模型液相扩散方程的新方法

谢奕展, 程夕明

A new method to solve electrolyte diffusion equations for single particle model of lithium-ion batteries

Xie Yi-Zhan, Cheng Xi-Ming
PDF
HTML
导出引用
计量
  • 文章访问数:  592
  • PDF下载量:  39
  • 被引次数: 0
出版历程
  • 收稿日期:  2021-09-01
  • 修回日期:  2021-10-20
  • 上网日期:  2022-02-13
  • 刊出日期:  2022-02-20

一种求解锂离子电池单粒子模型液相扩散方程的新方法

  • 北京理工大学机械与车辆学院, 电动车辆国家工程实验室, 北京 100081
  • 通信作者: 程夕明, cxm2004@bit.edu.cn
    基金项目: 国家重点研发计划(批准号: 2018YFB0106104)和国家自然科学基金(批准号: 51677006)资助的课题.

摘要: 电解液中的锂离子浓度表达是锂离子电池电化学模型求解的基本任务之一. 为了平衡单粒子模型的液相动态性能和计算效率, 假设反应仅发生在集电极和电解质界面上, 为此, 提出一种基于液相扩散方程无穷级数解析解的界面浓度求解新方法. 在恒流工况下, 利用数列单调收敛准则将解析解转化为一个收敛和函数. 在动态工况下, 将该解析解简化为输入与和函数的无限离散卷积. 利用和函数随时间单调衰减并收敛至零的特性对其进行截断, 从而得到有限离散卷积求解算法. 对比专业有限元分析软件, 该方法在恒流工况和动态工况下均能以较少的计算时间获得相当好的精度. 而且, 该方法仅有一个配置参数. 因此, 所提方法将有效减小应用于实时电池管理系统上的电化学模型计算负担.

English Abstract

    • 锂离子电池具有能量密度高、功率密度高、循环寿命长等特点[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液相扩散过程, 假设电化学反应仅发生在电解质与集电极界面处, 从而仅需求解界面处锂离子浓度. 为了高效、准确地获得恒流工况和动态工况下的界面锂离子浓度, 一种基于液相扩散方程解析解的有限离散卷积方法被提出.

    • 在改进单粒子模型中, 考虑液相的影响, 并假设正负极及隔膜的电导率、扩散系数等电化学性质相同, 且保持不变, 那么液相的质量守恒表达为

      $\begin{split} j(t,x) = \;&- D\nabla c(t,x) + \frac1{F}{{i(t,x){t_ + }}}, \\ & t \geqslant 0, ~~- L \leqslant x \leqslant L, \end{split} $

      $ \frac{{\partial c(t,x)}}{{\partial t}} + \frac{{\partial j(t,x)}}{{\partial x}} = 0, $

      $ i(t,x) = j(t,x)F, $

      式中, 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).

      电势守恒方程为

      $ i(t,x) {=} - \sigma \nabla \phi (t,x) {+} \frac{{2\sigma {R_g}T}}{F}(1 - {t_ + })\nabla \ln c(t,x), $

      式中, σ为液相电导率(S/m); ϕ(x, t)为液相厚度方向x处在t时刻的电势(V); Rg为气体常数; T为温度(K); 在改进单粒子模型中, 如图1液相模型所示, 假设固相与液相锂离子仅在液相边界处, 即x = –L, L处进入, 则边界条件为

      图  1  液相模型

      Figure 1.  The electrolyte model.

      $ j(t, - L) = j(t,L) = j(t), $

      式中, j(t)为正负极处固相与液相的边界锂离子交换通量密度(mol·m–2·s–1).

      假设初始时刻浓度处处分布均匀, 故初值条件为

      $ c(0,x) = {c_0}, $

      规定–L处电势为零点, 故另一边界条件为

      $ \phi (t, - L) = 0, $

      由于考虑液相的单粒子模型中固相与液相仅在边界处发生反应, 根据Bulter-Volmer方程, 液相仅需要求取液相两侧的锂离子浓度和电势, 也就是c(t, L), c(t, –L), ϕ (t, L)和ϕ (t, –L).

      联立(1)式、(2)式、(3)式、(5)式和(6)式, 得

      $ \left\{\begin{aligned} &\frac{\partial c(t,x)}{\partial t}-\frac{D}{1-{t}_+}\frac{{\partial }^{2}c(t,x)}{\partial {x}^{2}}=0, \\ & \qquad t > 0,~~ -L < x < L, \\ &c(0,x)={c}_{0}, \\ &-\frac{D}{1-{t}_+}\frac{{\rm{d}}c(t,-L)}{{\rm{d}}x}=j(t), \\ &-\frac{D}{1-{t}_+}\frac{{\rm{d}}c(t,L)}{{\rm{d}}x}=j(t) .\end{aligned} \right.$

      由(8)式可求解得到c(t, L)和c(t, –L).

      将(3)式代入(4)式, 并对(4)式进行积分, 得

      $\begin{split} &F\int_{ - L}^L {j(t,x){{\rm{d}}}x} = - \sigma [\phi (t,L) - \phi (t, - L)]\\ &+ \frac{{2\sigma {R_g}T}}{F}(1 - {t_ + })[\ln c(t,L) - \ln c(t, - L)]\;. \end{split}$

      由(1)式和(3)式联立, 得

      $ j(t,x) = - \frac{D}{{1 - {t_ + }}}\nabla c(t,x) . $

      将(7)式和(10)式代入(9)式中, 有

      $ \begin{split} &- \frac{{FD}}{{1 - {t_ + }}}[c(t,L) - c(t, - L)] = - \sigma [\phi (t,L)] \\ &+ \frac{{2\sigma {R_g}T}}{F}(1 - {t_ + })[\ln c(t,L) - \ln c(t, - L)] .\end{split} $

      由(11)式可知, ϕ (t, L)可由c(t, L)和c(t, –L)直接求得.

      综上, 改进单粒子模型液相部分仅需求解式(8), 从而得到液相两侧的浓度c(t, L)和c(t, –L).

    • 在(8)式中, 迁移系数t+能被整合进扩散系数D, 等效为一个新的扩散系数, 因此将t+视为0并不影响求解(8)式方法的一般性. 假设初始浓度为0, 那么在恒流工况下, 锂离子的液相扩散过程可由偏微分方程及其初始条件、边界条件描述:

      $ \left\{\begin{aligned} &\frac{\partial c(t,x)}{\partial t}-D\frac{{\partial }^{2}c(t,x)}{\partial {x}^{2}}=0 \\ & ~~ t > 0,~~~ -L < x < L, \\ &c(0,x)=0, \\ &-D\frac{{\rm{d}}c(t,-L)}{{\rm{d}}x}={j}_{\rm{const}}, \\ &-D\frac{{\rm{d}}c(t,L)}{{\rm{d}}x}={j}_{\rm{const}}, \end{aligned} \right.$

      式中, jconst为液相两侧界面通量密度(mol·m–2·s–1).

      对(12)式进行拉普拉斯变换, 得

      $ \left\{\begin{aligned} &D\frac{{\partial }^{2}C(x)}{\partial {x}^{2}}-sC(x)=0, \\ &-D\frac{\partial C(-L)}{\partial x}=\frac{{j}_{\rm{const}}}{s}, \\ &-D\frac{\partial C(L)}{\partial x}=\frac{{j}_{\rm{const}}}{s}, \end{aligned} \right.$

      式中, C(x)是c(t, x)的拉普拉斯变换. 求解(13)式, 得

      $ \begin{split} &C(x) = {j_{{\rm{const}}}}\frac{1}{s}\frac{1}{{\sqrt {sD} }}\\ &\times\frac{{\cosh [\sqrt {s/D} (x - L)] - \cosh [\sqrt {s/D} (x + L)]}}{{\sinh (2\sqrt {s/D} L)}} .\end{split} $

      则液相两侧浓度的拉普拉斯变换为

      $ C(-L)={j}_{\rm{const}}\frac{1}{s}\frac{1}{\sqrt{sD}}\frac{{\rm{cosh}}[\sqrt{s/D}R]-1}{{\rm{sinh}}(\sqrt{s/D}R)}, $

      $ C(L) = - {j_{{\rm{const}}}}\frac{1}{s}\frac{1}{{\sqrt {sD} }}\frac{{\cosh (\sqrt {s/D} R) - 1}}{{\sinh (\sqrt {s/D} R)}}\;. $

      显然, C(L)和C(–L)互为相反数, 因此在下文内容中将仅考虑C(L). 对(16)式进行拉普拉斯逆变换, 得L侧的浓度为

      $\begin{split} & {c_L}(t) = - {j_{{\rm{const}}}}\frac{R}{{2D}} + {j_{{\rm{const}}}}\frac{{4R}}{D}\mathop \sum \limits_{n = 1}^\infty \bigg\{ \frac{1}{{{{(2n - 1)}^2}{{{\rm{\pi}}}^2}}} \exp [ - {(2n - 1)^2}{{{\rm{\pi}}}^2}\Big(\dfrac{D}{R^2}\Big)t]\bigg\} \;. \end{split} $

      由于(17)式是一无穷级数, 需要进行简化以方便求解. 令收敛和函数:

      $ f(t) = \mathop {\lim }\limits_{N \to \infty } \mathop \sum \limits_{n = 1}^N \left\{ \frac{1}{{{{(2n - 1)}^2}{{{\rm{\pi}}}^2}}}\exp [ - {(2n - 1)^2}{{{\rm{\pi}}}^2}\Big(\dfrac{D}{R^2}\Big)t]\right\} \;. $

      并令

      $ {f_N}(t) = \mathop \sum \limits_{n = 1}^N \left\{ \frac{1}{{{{(2n - 1)}^2}{{{\rm{\pi}}}^2}}}\exp [ - {(2n - 1)^2}{{{\rm{\pi}}}^2}\Big(\dfrac{D}{R^2}\Big)t]\right\} \;. $

      由于

      $ \mathop {\lim }\limits_{N \to \infty } {f_N}(0) = \mathop {\lim }\limits_{N \to \infty } \mathop \sum \limits_{n = 1}^N \frac{1}{{{{(2n - 1)}^2}{{{\rm{\pi}}}^2}}} = \frac{1}{8}\;. $

      且{fN(0)}是一单调递增数列, 所以

      $ {f_N}(0) \leqslant \frac{1}{8},~~\forall~ N . $

      又因为fN(t) > 0且随t单调衰减, 所以对任意t > 0, 数列{fN(t)}均单调递增且有上下界. 从而, {fN(t)}存在极限, 即(18)式的f(t)存在且是一有界函数.

      因此, (17)式可重写为

      $ {c_L}(t) = {j_{{\rm{const}}}}\frac{{4R}}{D}\left[ {f(t) - \frac{1}{8}} \right]\;. $

      $ g(t) = f(t) - \frac{1}{8} . $

      那么, (22)式可表达为

      $ {c_L}(t) = {j_{{\rm{const}}}}\frac{{4R}}{D}g(t) . $

    • 在动态工况下, (12)式改写为

      $ \left\{\begin{aligned} &\frac{\partial c(t,x)}{\partial t}-D\frac{{\partial }^{2}c(t,x)}{\partial {x}^{2}}=0, \\ & \qquad t > 0,~~ -L < x < L, \\ &c(0,x)=0, \\ &-D\frac{{\rm{d}}c(t,-L)}{{\rm{d}}x}=j(t), \\ &-D\frac{{\rm{d}}c(t,L)}{{\rm{d}}x}=-j(t), \end{aligned} \right.$

      式中 j(t)为液相两侧端点时变通量密度(mol·m–2·s–1). 由(24)式得系统的单位脉冲响应为

      $ {\delta _L}(t) = \frac{{4R}}{D}g'(t) . $

      则液相L侧的浓度为

      $ {c_L}(t) = \frac{{4R}}{D}\int_0^t {j(t - \tau )g'(\tau ){{\rm{d}}}\tau } . $

      对(27)式进行分部积分, 得

      $ {c_L}(t) = \frac{{4R}}{D}\left[j(0)g(t) - \int_0^t {j'(t - \tau )g(\tau ){{\rm{d}}}\tau } \right] . $

      $ y(t) = \int_0^t {j'(t - \tau )g(\tau ){{\rm{d}}}\tau } . $

      则(28)式可重写为

      $ {c_L}(t) = \frac{{4R}}{D}[j(0)g(t) - y(t)] . $

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

    • (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, 则有

      $ {j_n} = {j_{n - 1}} + {k_n}{{\rm{d}}}t, $

      $ j({t_{n - 1}} + \tau ) = {j_{n - 1}} + {k_n}\tau, $

      式中, knj(t)在区间[tn-1, tn]上的斜率. 从而, (29)式可改写为

      $ \begin{split} y(t) =\;& \int_0^t {j'(t - \tau )g(\tau ){{\rm{d}}}\tau } \\ =\;& - \sum\limits_{n = 1}^N {[{k_{N + 1 - n}}\int_{{t_{n - 1}}}^{{t_n}} {g(\tau ){\rm{d}}\tau } ]} . \end{split} $

      $ \begin{split} \;& \int_0^t {g(\tau ){\rm{d}}\tau } = \int_0^t {\left[ {\left( {\mathop {\lim }\limits_{N \to \infty } \sum\limits_{n = 1}^N {\left\{ {\frac{{\exp [ - {{(2n - 1)}^2}{\pi ^2}(D/{R^2})\tau ]}}{{{{(2n - 1)}^2}{\pi ^2}}}} \right\}} - \frac{1}{8}} \right){\rm{d}}\tau } \right]} \\ = \;&- \frac{1}{8}t + {\kern 1pt} \mathop {\lim }\limits_{N \to \infty } \mathop \sum \limits_{n = 1}^N \frac{1}{{{{(2n - 1)}^2}{{{\rm{\pi}}}^2}}}\int_0^t {\exp [ - {{(2n - 1)}^2}{{{\rm{\pi}}}^2}(D/{R^2})\tau ]{{\rm{d}}}\tau } \\ = \;&- \frac{1}{8}t + \frac{{{R^2}}}{D}\mathop {\lim }\limits_{N \to \infty } \mathop \sum \limits_{n = 1}^N \frac{1}{{{{(2n - 1)}^4}{{{\rm{\pi}}}^4}}} - \frac{{{R^2}}}{D}\mathop {\lim }\limits_{N \to \infty } \mathop \sum \limits_{n = 1}^N \frac{{\exp [ - {{(2n - 1)}^2}{{{\rm{\pi}}}^2}(D/{R^2})\tau ]}}{{{{(2n - 1)}^4}{{{\rm{\pi}}}^4}}} \end{split} . $

      记另一收敛和函数:

      $ h(t) = \mathop {\lim }\limits_{N \to \infty } \mathop \sum \limits_{n = 1}^N \frac{{\exp [ - {{(2n - 1)}^2}{{{\rm{\pi}}}^2}(D/{R^2})\tau ]}}{{{{(2n - 1)}^4}{{{\rm{\pi}}}^4}}} . $

      f(t)的分析类似, h(t)极限存在且单调有界. 简记hn = h(tn), 有

      $ \int_{{t_{n - 1}}}^{{t_n}} {g(\tau ){{\rm{d}}}\tau } = - \frac{1}{8}{{\rm{d}}}t - \frac{{{R^2}}}{D}({h_n} - {h_{n - 1}}) . $

      将(36)式代入(33)式中, 有

      $ \begin{split} y(t) =\;& \sum\limits_{n = 1}^N {\left\{ {{k_{N + 1 - n}}\left[ {\frac{1}{8}{\rm{d}}t + \frac{{{R^2}}}{D}({h_n} - {h_{n - 1}})} \right]} \right\}} \\ =\;& \frac{{{R^2}}}{D}\sum\limits_{n = 1}^N {[{k_{N + 1 - n}}({h_n} - {h_{n - 1}})} ] \\ &+ \frac{1}{8}{{\rm{d}}}t\sum\limits_{n = 1}^N {{k_{N + 1 - n}}}.\\[-15pt] \end{split} $

      将(37)式代入(30)式中, 并简记g(tn) = gn, 整理得动态工况下浓度的离散卷积形式解析解为

      $\begin{split} {c_L}({t_N}) = \;&\frac{{4R}}{D}\left\{ {{j_0}{g_n} - \frac{{{R^2}}}{D}\sum\limits_{n = 1}^N {[{k_{N + 1 - n}}({h_n} - {h_{n - 1}})} ] - \frac{1}{8}{\rm{d}}t\sum\limits_{n = 1}^N {{k_{N + 1 - n}}} } \right\}\\ =\;& \frac{{4R}}{D}{j_0}{g_N} - \frac{R}{{2D}}({j_N} - {j_0})\; - 4\frac{{{R^3}}}{{{D^2}}}\sum\limits_{n = 1}^N {[{k_{N + 1 - n}}({h_n} - {h_{n - 1}})} ]. \end{split} $

    • 显然(38)式的计算量随着t的增大而线性增长, 使得其在t很大时不具备可行性, 需要对其进行简化. 注意到随h(t)随t单调衰减, 且有

      $ \mathop {\lim }\limits_{t \to \infty } h(t) = 0 . $

      根据表1参数绘制得h(t)的函数图像如图2所示. 显然, 当t大于某个值, 称为截断时长T, 且假设此时N0 = T/dt, 可认为nN0时, hnhn–1 ≈ 0, 并考虑初始浓度, 由(38)式简化得浓度的有限离散卷积算法为:

      参数单位
      扩散系数D7.5 × 10–11m2/s
      液相总厚度R3.35 × 10–4m
      初始浓度c02000mol/m3
      1 C放电时液相两侧界面通量 j1C1.74 × 10–4mol·m–2·s–1
      *ComSol

      表 1  模型参数*

      Table 1.  Model Parameters*.

      图  2  h(t)函数曲线

      Figure 2.  The curve of h(t).

      $ \begin{split} &{c_L}(N{{\rm{d}}}t) = {c_0} + \frac{{4R}}{D}{j_0}{g_N} - \frac{R}{{2D}}({j_N} - {j_0}) \\ & - \left\{ {\begin{aligned} &{4\frac{{{R^3}}}{{{D^2}}}\sum\limits_{n = 1}^N {[{k_{N + 1 - n}}({h_n} - {h_{n - 1}})]} }, &{N < {N_0}} \\ &{4\frac{{{R^3}}}{{{D^2}}}\sum\limits_{n = 1}^{{N_0}} {[{k_{N + 1 - n}}({h_n} - {h_{n - 1}})]} }, &{N \geqslant {N_0}} . \end{aligned}} \right. \end{split} $

    • 由于放电和充电工况的液相两侧表面通量的方向相反, 因此以放电和充放电混合工况分析求解方法的性能. 一般地, 纯电动汽车动力电池系统的最大放电倍率不超过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), 为最大绝对误差与初始浓度之比.

      $ {\rm{RRMSE}} = \frac{1}{{{c_0}}}\sqrt {\frac{1}{m}\mathop \sum \limits_{i = 1}^m {{({u_i} - {v_i})}^2}}, $

      $ {\rm{RMAE}} = \frac{1}{{{c_0}}}{\max _i}|{u_i} - {v_i}|, $

      式中, c0表1中的初始浓度; m为样本数目; ui为所提方法计算值; vi为参考解.

    • 在恒流工况下, 电极两侧浓度表达式为(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误差

      Figure 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)函数以及相关系统参数, 计算简便.

    • 在动态工况下, 电极两侧锂离子浓度由(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) 误差

      Figure 4.  Simulation results for the DST profile: (a) Profile; (b) concentration; (c) error.

      图  5  NEDC工况仿真结果 (a) 放电倍率; (b) 浓度; (c) 误差

      Figure 5.  Simulation results for the NEDC profile: (a) Profile; (b) concentration; (c) error.

      图  6  WLTC工况仿真结果 (a) 放电倍率; (b) 浓度; (c) 误差

      Figure 6.  Simulation results for the WLTC 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仿真结果高度一致. 出现上述误差现象是由图1h(t)在不同截断时长截断所产生. 显然, N0越小, 截断时长T越小, h(T)越小, 导致(38)式截断为(40)式所忽略的项:

      $ \sum\limits_{n = {N_0}}^N {[{k_{N + 1 - n}}({h_n} - {h_{n - 1}})]} . $

      越大, 使截断误差增大. 而当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%.

      参数RRMSE(%)RMAE(%)
      AnaN0(200)DST0.99552.8533
      AnaN0(∞)DST0.08800.6093
      AnaN0(200)NEDC1.87874.882
      AnaN0(600)NEDC0.09910.3123
      AnaN0(1000)NEDC0.03730.1635
      AnaN0(∞)NEDC0.03830.1779
      AnaN0(200)WLTC2.67676.0554
      AnaN0(600)WLTC0.14050.3221
      AnaN0(1000)WLTC0.04720.1425
      AnaN0(∞)WLTC0.04630.1565

      表 2  动态工况误差

      Table 2.  Errors under dynamic profiles.

    • 以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增大时, 所需计算的$ k_{N + 1 - n}(h_n - h_{n-1}) $项增多, 导致求解时间有所增长, 但相对求解时间都不超过1/1000. 在N0 = 600时, 卷积算法的结果AnaN0(600)WLTC与Comsol结果eT(1 × 10–3)WLTC相比, 二者求解精度接近, 但所提卷积算法的求解时间仅为Comsol的1/500.

      参数RRMSE(%)RT(%)
      eT(1 × 10–3)WLTC0.121835.14
      AnaN0(200)WLTC2.67670.06166
      AnaN0(600)WLTC0.14050.06765
      AnaN0(1000)WLTC0.04720.09490
      AnaN0($\infty $)WLTC0.04630.09573

      表 3  WLTC工况下卷积算法与Comsol比较

      Table 3.  The comparison between the convolution algorithm and COMSOL under the WLTC.

    • 本文提出了一种新的求解锂离子电池电化学模型液相扩散过程液相两侧表面浓度的离散卷积算法. 该算法求解简单, 参数较少, 仅需配置1个参数, 模型的求解误差随着截断时长T的增加而减小.

      该离散卷积算法是基于扩散方程的解析解而提出, 具有收敛快速而精度高的优点. 在DST, NEDC和WLTC工况下, 相比于容许误差eT = 1 × 10–5的Comsol模型求解结果, 所提算法取截断误差最大时的计算误差也仅为3%, 而取截断误差最小时的计算误差仅为0.04%左右, 具有相当高的计算精度. 在WLTC工况激励下, 该卷积算法能够以有限元方法1/500的求解时间取得与后者相近的精度.

      在数据存储成本方面, 该离散卷积算法仅需存储工况数据、收敛和函数g(t)和h(t)的各N0个样本值, 数据存储量少. 因此, 所提方法将有助于电化学模型在电池管理系统上的实时应用. 对于不同的电池液相参数, 特别是时变参数, 算法的适用性需进一步研究. 在下一阶段工作中, 将聚焦于该离散卷积算法在锂离子电池电化学降阶模型上的集成.

参考文献 (21)

目录

    /

    返回文章
    返回