一种全工况多尺度动力锂电池电化学耦合建模方法

1.本发明属于锂离子电池技术领域，尤其涉及一种全工况多尺度动力锂电池电化学耦合建模方法。


背景技术：

2.随着电动汽车产业的迅速崛起，动力锂电池等高能量密度电源的需求也日益增加。
3.动力锂电池仿真模型主要有等效电路模型和电化学模型，等效电路模型利用简单的电器元件来合理搭建电路网络，以此模拟电池的一些基础结构和物理化学过程的外在表现，难以反应电池内部物理、化学过程的实际运行情况。锂电池是一个多维度宽物理层的系统，电化学模型便从电池外部和内部物理特性出发构造模拟电池，充分考虑了电池内部的实际电化学反应过程，包括固液相锂离子的扩散和迁移过程，固液相电势和电荷守恒问题等，完整描述电池内反应机理，与其他模型相比具有更高的精度。
4.经典的p2d模型精度较高，但其耦合复杂的非线性偏微分方程，相对于等效电路模型具有较高的复杂度，求解困难、运算效率低，这就限制了其在控制系统中的应用。单粒子模型相对于p2d模型计算效率有所提高，但由于忽略液相动力学行为，在高倍率放电情况下误差较大，精度无法保证，使其无法直接应用于工程实际中，因此电化学模型的进一步降阶简化在电池管理中十分重要。


技术实现要素：

5.本发明所要解决的技术问题是提供一种全工况多尺度动力锂电池电化学耦合建模方法来解决现有技术中经典p2d模型过于复杂，计算量过大，而单粒子模型中又缺少对液相动力学行为的考虑，模型精度低，使得这两种电化学模型都无法直接应用于电池管理系统中的问题。为了平衡p2d模型的液相动态性能、计算效率和普适性，本发明提出一种在全工况充放电条件下，基于大系统理论简化后得到的简化固相扩散方程和基于有限离散卷积方法求解得到的简化液相扩散方程相耦合，提高模型运算速度的同时增强模型的适用性。锂离子在单粒子颗粒物内的扩散过程是沿球形颗粒径向维度，同时还需要考虑时间维度的信息，因此该发明称为全工况多尺度动力锂电池电化学耦合模型。
6.为解决上述技术问题，本发明的技术方案为：提供一种全工况多尺度动力锂电池电化学耦合建模方法，其创新点在于：具体包括以下步骤：
7.步骤1：设固相锂离子浓度为cs、锂离子孔壁流量密度为j，设固相锂离子浓度cs和锂离子孔壁流量密度j在沿电池负极到正极的方向上的任一时刻是常量；对固相浓度cs的扩散方程进行非均匀离散化处理，离散点数为mr，得到抛物线参数方程，来创建若干个分层半径rk，若干个分层半径rk再结合固相浓度扩散方程得到一个输入为电流i，输出是颗粒表面锂离子浓度css的状态空间方程；
8.步骤2：将步骤1中得到的状态空间方程中的系数矩阵划分为块矩阵，分解出显性
模块和非显性模块以选取特征值，确定n个主导状态变量，并按照顺序变换矩阵得到1
×
n阶的子模型矩阵；计算步骤1中系数矩阵的特征值与特征向量，并从特征向量中构造矩阵λi(i＝1,2,
…mr-1),vi＝(1,2,
…mr-1)；
9.步骤3：计算主要标准并选择步骤2中n个具有主要标准的主要特征值，基于λi和vi得到的新阶变换矩阵v，g，划分为两部分块矩阵；计算误差矩阵e，根据计算得到的数据，选择显性特征值，经过大系统理论简化后构造出简化状态空间方程，即简化固相模型；
10.步骤4：根据bulter-volmer方程，推导得到液相两侧的锂离子浓度与时间的函数关系式；
11.步骤5：在恒流工况下求解步骤4中的锂离子浓度扩散方程；
12.步骤6：在动态工况下求解步骤4中的锂离子浓度扩散方程，通过离散化处理得到简化液相模型；
13.步骤7：根据工况要求选取步骤5或步骤6中求解得到的正负极液相锂离子浓度，在整个电池的x轴方向上对液相电势微分方程进行积分，可得到正负极液相电势的差值；
14.步骤8：取步骤3中的固相锂离子浓度，根据工况要求选取步骤5或步骤6中的液相锂离子浓度，求解bulter-volmer程，得到正负极过电势表达式；
15.步骤9：取步骤3中计算得到的颗粒表面锂离子浓度，求解正负极开路电压之差，最后得到简化的端电压计算公式，完成全工况多尺度动力锂电池电化学耦合建模工作。
16.进一步的，所述的步骤1中假设的锂离子孔壁流量密度j为：
[0017][0018]
其中，i是电池工作电流密度，f是法拉第常数，a是电极单位体积有效反应面积，l是沿x轴的不同区域电池厚度， 和-分别代表正负电极；
[0019]
所述的步骤1中，固相浓度cs的扩散方程为：
[0020][0021]
对固相浓度cs的扩散方程进行非均匀离散化处理，得到的抛物线参数方程为：
[0022]
y＝[-(k-mr)2 mr2]
×rs
/mr2[0023]
若干个分层半径rk结合固相浓度扩散方程得到的状态空间方程为：
[0024][0025]
其中，a2、b2、c2、d2都是系数矩阵，a2是非奇异矩阵，具体为：
[0041]
最后简化的子模型可以表示为
[0042]
进一步的，所述的步骤3中设主要标准为mk和sk，则主要标准mk和sk的计算方法为：
[0043]
对z
′
＝∧z gu进行拉普拉斯变换，得到基本状态变量x的传递函数，并分析系统的动态响应，当输入信号为输出信号为x1(s),x
1j
(s)表示状态变量的相应响应，时域响应为x
1j
(t)：
[0044][0045][0046][0047]
设cj,k表示由引起的响应幅度：
[0048]cj,k
＝v
j,kgk
/λk[0049]
输入与每个基本状态变量之间一一对应，所以定义如下mk和sk来分析每一个特征值的显著性，并据此选取显性特征值:
[0050][0051][0052]
其中，上述两个公式即为主要标准，所述的步骤3中，对z
′
＝∧z gu进行拉普拉斯变换,当输入信号为阶跃u＝1/s时，显性和非显性模式的时域响应表示为z1(t)，z2(t)；εi与误差矩阵e的行向量相关，基于此提出一个成本函数：
[0053][0054]
则目标函数表示为:
[0055]
[0056][0057]
ηi是拉格朗日乘数(i＝1，2，
…
，n-m)；
[0058]
步骤3中所述的误差矩阵e为：
[0059][0060][0061]
进一步的，步骤4中所述的液相两侧的锂离子浓度与时间的函数关系式为：
[0062][0063]
p2d模型中固相与液相仅在边界处发生反应，根据bulter-volmer方程可知，求解端电压时，液相中只需要求解两侧的锂离子浓度与电势，de是液相扩散系数，t

是离子转移数。
[0064]
进一步的，所述的步骤5中恒流工况下液相浓度扩散方程求解过程如下：
[0065]
假设初始浓度为0，将t

整合到新的扩散系数d中，对步骤4中所述的液相两侧的锂离子浓度与时间的函数关系式进行拉普拉斯变换和拉普拉斯逆变换得到两侧的锂离子浓度：
[0066][0067]
表示液相两侧的锂离子浓度，经计算可知，f(t)是单调递增且有上下界的一个收敛和函数：
[0068][0069]
[0070]jconst
是液相两侧界面通量密度(mol/m2/s)，r为液相区域总厚度，其中和互为相反数。
[0071]
进一步的，步骤6中所述的动态工况下液相浓度扩散方程求解过程如下：
[0072]
从拉普拉斯变换和拉普拉斯逆变换后得到的锂离子浓度与时间的关系式中，得到系统单位脉冲响应，接着对液相扩散方程进行分步积分，得到：
[0073][0074][0075]
液相两侧的锂离子浓度与时间的函数关系式进行拉氏变换，将区间[0,t]等分为n个区间，第n个区间为[t
n-1
,tn]，tn＝ndt,每个区间大小为dt＝t/n，以1秒为一个采样点，假设节点之间的值由线性插值获得，简记j(tn)＝jn,
[0076]jn
＝j
n-1
kndt
[0077]
j(t
n-1
τ)＝j
n-1
knτ
[0078]
得到：
[0079][0080]kn
是j(t)在区间[t
n-1
,tn]上的斜率，其中令收敛和函数为h(t)，经计算得知h(t)的极限存在且是单调有界函数，简记hn＝h(tn):
[0081][0082]
则有
[0083][0084]
简记g(tn)＝gn，整理上述公式，得到动态工况下液相锂离子浓度扩散方程的离散卷积解为：
[0085][0086]
步骤6中得到的离散卷积解需要进行截断简化处理以提高运算效率：
[0087]
规定t大于某个值时存在截断时长t并且n0＝t/dt；在考虑有液相动力学影响的初始液相锂离子浓度时，n≥n0则有h
n-h
n-1
≈0，故简化后的液相离散卷积算法为：
[0088][0089]
进一步的，步骤7中所述的液相电势之差的求解过程：
[0090]
令且在x轴上为常量，
[0091][0092]
对上式在电池的x轴方向上进行积分，得到正负极电势之差为：
[0093][0094]keff
是电解液有效电导率，r是气体常数，t是参考温度，f
c/a
是液相平均摩尔活度系数，sep代表电极隔膜。
[0095]
进一步的，所述的步骤8中的过电势表达式为：
[0096][0097]
其中
[0098]
i0±
＝k[c
±
lcss
(c
s,max-c
ss
)]
α
[0099]
α是电极反应转化系数，i
0
是正负极集流板电流密度，a是正负极单位体积有效反应面积，k是电解液电导率。
[0100]
进一步的，步骤9中所述的端电压计算表达式：
[0101][0102]
和现有技术相比，本发明的有益效果为：
[0103]
(1)本发明首次通过大系统理论和有限离散卷积两种方法，分别对经典p2d模型的固液相扩散方程降阶简化之后进行耦合，相对于经典p2d模型而言，本发明模型在保证精度
的同时大大减小了模型的复杂度。
[0104]
(2)本发明考虑了液相动力学行为的影响，在继承单粒子模型高运算速率的优点的同时，大幅提高了模型的精度，且模型精度不受充放电倍率的影响。
[0105]
(3)本发明通过假设锂离子孔壁流量密度j在沿电池负极到正极的方向上的任一时刻是常量，对液相电势微分方程进行全积分，得到简化的电势表达式，进一步降低整体模型的复杂度。
[0106]
(4)本发明在恒流和动态全工况下融合简化的固液相简化模型，首次得到一种适应性强且运算时间短的高精度耦合模型，可应用于多工况电池管理系统中。
附图说明
[0107]
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。
[0108]
图1是本发明非均匀离散差分的实施图。
[0109]
图2是液相扩散方程的液相模型。
[0110]
图3是本发明实施的流程示意图。
具体实施方式
[0111]
为了说明本发明所述的技术方案，下面结合附图和公式进一步说明本发明。
[0112]
本发明提供一种全工况多尺度动力锂电池电化学耦合建模方法，其创新点在于：具体包括以下步骤：
[0113]
步骤1：设固相锂离子浓度为cs、锂离子孔壁流量密度为j，设固相锂离子浓度cs和锂离子孔壁流量密度j在沿电池负极到正极的方向上的任一时刻是常量，假设的锂离子孔壁流量密度j为：
[0114][0115]
其中，i是电池工作电流密度，f是法拉第常数，a是电极单位体积有效反应面积，l是沿x轴的不同区域电池厚度， 和-分别代表正负电极；
[0116]
固相浓度cs的扩散方程为：
[0117][0118]
对固相浓度cs的扩散方程进行非均匀离散化处理，离散点数取mr，得到抛物线参数方程来创建若干个分层半径rk，得到mr关于y的抛物线参数方程，如图1所示，横坐标是离散点数，纵坐标是抛物线参数方程的值，锂离子颗粒直径等于方程的最大值，非均匀离散得到的颗粒半径与抛物线参数方程的值一一对应：
[0119]
y＝[-(k-mr)2 mr2]
×rs
/mr2[0120]
若干个分层半径rk结合固相浓度扩散方程得到一个输入为电流i，输出为颗粒表面锂离子浓度css的状态空间方程：
[0121][0122]
其中，a2、b2、c2、d2都是系数矩阵，a2是非奇异矩阵，具体为：
[0123][0124][0125][0126][0127]
αk＝2ds/(δr
k-1
δrk)
[0128]
δrk＝r
k 1-rk[0129]
其中k＝2
…mr-2
，ds是固相浓度扩散系数；
[0130]
步骤2：将步骤1中得到的状态空间方程中的系数矩阵划分为块矩阵，分解出显性模块和非显性模块以选取特征值，确定n个主导状态变量，并按照顺序变换矩阵得到1
×
n阶的子模型矩阵；计算步骤1中系数矩阵的特征值与特征向量，并从特征向量中构造矩阵λi(i＝1,2,
…mr-1),vi＝(1,2,
…mr-1)；其中，块矩阵推导过程如下：z
′
＝∧z gu
[0131]
∧＝∨-1
[0132]
[0133][0134]z′1＝∧1z1 g1u 显性模块
[0135]z′2＝∧2z2 g2u 非显性模块
[0136]
x＝fz1[0137]
最后简化的子模型可以表示为
[0138]
步骤3：设主要标准为mk和sk，计算主要标准并选择步骤2中n个具有主要标准的主要特征值，主要标准mk和sk的计算方法为：对z
′
＝∧z gu进行拉普拉斯变换，得到基本状态变量x的传递函数，并分析系统的动态响应，当输入信号为输出信号为x1(s),x
1j
(s)表示状态变量的相应响应，时域响应为x
1j
(t)：
[0139][0140][0141][0142]
设cj,k表示由引起的响应幅度：
[0143]cj,k
＝v
j,kgk
/λk[0144]
输入与每个基本状态变量之间一一对应，所以定义如下mk和sk来分析每一个特征值的显著性，并据此选取显性特征值:
[0145][0146][0147]
其中，上述两个公式即为主要标准；
[0148]
基于λi和vi得到的新阶变换矩阵v，g，划分为两部分块矩阵：对z
′
＝∧z gu进行拉普拉斯变换,当输入信号为阶跃u＝1/s时，显性和非显性模式的时域响应表示为z1(t)，z2(t)；εi与误差矩阵e的行向量相关，基于此提出一个成本函数：
[0149][0150]
即
[0151][0152][0153][0154]
则目标函数表示为:
[0155][0156][0157]
ηi是拉格朗日乘数(i＝1，2，
…
，n-m)；
[0158]
计算误差矩阵e，步骤3中所述的误差矩阵e为：
[0159][0160][0161]
根据计算得到的数据，选择显性特征值，经过上述步骤简化后构造出简化状态空间方程，即简化固相模型，因为构造出简化状态空间方程时，取的n值不确定，显性特征值也就随机的，所以简化状态空间方程不固定。
[0162]
步骤4：图2所示的液相模型中，j(t)是正负极处固相与液相的边界锂离子交换通量密度，the electrolyte部分是电解液，根据bulter-volmer方程和图2所示的液相模型，推导得到液相两侧的锂离子浓度与时间的函数关系式为：
[0163][0164]
p2d模型中固相与液相仅在边界处发生反应，根据bulter-volmer方程可知，求解端电压时，液相中只需要求解两侧的锂离子浓度与电势，c
l
(t,x)是液相锂离子浓度，de是液相扩散系数，t

是离子转移数；
[0165]
步骤5：在恒流工况下求解步骤4中的锂离子浓度扩散方程，求解过程如下：
[0166]
假设初始浓度为0，将t

整合到新的扩散系数d中，对步骤4中所述的液相两侧的锂离子浓度与时间的函数关系式进行拉普拉斯变换和拉普拉斯逆变换得到两侧的锂离子浓度：
[0167][0168]
表示液相两侧的锂离子浓度，其中和互为相反数，j
const
是液相两侧界面通量密度，r为液相区域总厚度，经计算可知，f(t)是单调递增且有上下界的一个收敛和函数：
[0169][0170][0171]
步骤6：在动态工况下求解步骤4中的锂离子浓度扩散方程，通过离散化处理得到简化液相模型，求解过程如下：
[0172]
从拉普拉斯变换和拉普拉斯逆变换后得到的锂离子浓度与时间的关系式中，得到系统单位脉冲响应，接着对液相扩散方程进行分步积分，得到：
[0173][0174][0175]
液相两侧的锂离子浓度与时间的函数关系式进行拉氏变换，将区间[0,t]等分为n个区间，第n个区间为[t
n-1
,tn]，tn＝ndt,每个区间大小为dt＝t/n，以1秒为一个采样点，假
设节点之间的值由线性插值获得，简记j(tn)＝jn,
[0176]jn
＝j
n-1
kndt
[0177]
j(t
n-1
τ)＝j
n-1
knτ
[0178]
得到：
[0179][0180]kn
是j(t)在区间[t
n-1
,tn]上的斜率，其中令收敛和函数为h(t)，经计算得知h(t)的极限存在且是单调有界函数，简记hn＝h(tn):
[0181][0182]
则有
[0183][0184]
简记g(tn)＝gn，整理上述公式，得到动态工况下液相锂离子浓度扩散方程的离散卷积解为：
[0185][0186]
步骤6中得到的离散卷积解需要进行截断简化处理以提高运算效率：
[0187]
规定t大于某个值时存在截断时长t并且n0＝t/dt；在考虑有液相动力学影响的初始液相锂离子浓度时，n≥n0则有h
n-h
n-1
≈0，故简化后的液相离散卷积算法为：
[0188][0189]
步骤7：取步骤6中求解得到的正负极液相锂离子浓度，在整个电池的x轴方向上对液相电势微分方程进行积分，可得到正负极液相电势的差值，液相电势之差的求解过程：
[0190]
令且在x轴上为常量，
[0191][0192]
对上式在电池的x轴方向上进行积分，得到正负极电势之差为：
[0193][0194]
l是电极厚度，k
eff
是电解液有效电导率，是电解液电导率，r是气体常数，t是参考温度，f
c/a
是液相平均摩尔活度系数，sep代表电极隔膜；
[0195]
步骤8：取步骤3～步骤6中计算得到的固液相锂离子浓度，求解bulter-volmer程，得到正负极过电势表达式为：
[0196][0197]
其中
[0198]
i0±
＝k[c
±
lcss
(c
s,max-c
ss
)]
α
[0199]
α是电极反应转化系数，i
0
是正负极集流板电流密度，a是正负极单位体积有效反应面积，k是电解液电导率；
[0200]
步骤9：取步骤3～步骤4中计算得到的颗粒表面锂离子浓度，求解正负极开路电压之差，最后得到简化的端电压计算公式：
[0201][0202]
图3是发明实施的流程示意图，同时也是全工况多尺度动力锂电池电化学耦合模型，从(1)(2)(3)三个方向对经典p2d模型使用上述发明案例中的相应方法进行简化降阶，最后进行耦合得到全工况多尺度动力锂电池电化学耦合模型。