Matlab&Comsol联合优化磷酸铁锂电池电化学产热模型参数

详细信息

作者优秀: 优秀教师/意见领袖/博士学历/特邀专家/独家讲师
平台推荐: 内容稀缺

16天前浏览2066

导读：一个好的仿真必定是简洁精确的理论框架和准确的实验数据相结合的。不能一味追求模型的复杂程度而忽视其应用效果，或者缺乏足够的数据支撑。因此建模须遵循规范的建模流程，做到环环相扣。本文简单演示对于磷酸铁锂电池，通过使用Matlab算法处理实验数据，并传递至Comsol进行仿真计算的案例，该案例后续也会以视频教程方式放进课程。推荐大家关注我的视频课程《电池系统管理Matlab&Simulink&Comsol联合仿真入门》。

一、方壳电池热安全问题

20Ah的磷酸铁锂电池由于体积较大，内部温度分布可能明显不均匀，对于其长时安全使用造成安全隐患。因此，建立一个准确的模型计算电池内部温度分布可以有效提升电池的安全寿命。本文在实验测量热学数据并使用Matlab进行参数识别后，使用Comsol计算其电化学模型，计算内部温度分布等参数。输入基于数据识别后的参数后，模型可以确地预测了动态终端电压，以及电池外部的最小、最大温度和表面平均温度。通过算法提取出的电化学和热学参数与文献参数较为吻合，表明可以通过算法拟合外部测量数据来推断材料的内部特性。

二、实验设计

电池的基础特性取决于其SOC状态。因此需要设计一种固定充电周期的充放电方案，在不引起热力学状态变化的情况下研究速率和温度的关系。在SOC=30%和50%的情况下，使用不同幅度与周期的方波循环，在该SOC区间电池电压曲线在一个平台，内部参数相对稳定。确定由输入信号引起的SOC摆动后，初始电池的SOC进行偏移，使充电期间的平均状态以所选择的SOC为均值。例如，对于4C倍率，100s为周期的方波，充电部分使电池SOC提升5.56%，为了使循环期间的平均SOC=30%，电池的初始SOC需要维持在30-5.56/2=27.22%。电池被专门设计的夹具固定，该夹具可以保持电池热学性能的稳定。

电池进行交流激励时，表面温度分布是用红外相机测量的。实验是在一个热隔离的房间里进行，稳定的环境温度由外部空调控制。

方波激流的电压相应和温度测量曲线如下图所示

图中曲线显示电池表面的最小值、最大值和表面平均温度的变化。最大温度波动ΔT也被视为特征值并用于参数化，用于量化温度分布的非均匀性。

三、Comsol模型设置及结果

电荷守恒方程：

欧姆定律

产热方程

模型几何

温度变化

温度曲线

四、Matlab模型设置

为了找到热参数和电化学参数的最优值，通过拟合周期稳态数据（从~2000秒开始）进行非线性最小二乘优化，包括4c100秒方波的电池电压、最小值、最大值和平均温度。采用MATLAB最小二乘法，取目标函数为：

其中，上标和表示模型和实验结果，和分别为参数的下界和上界；表示模拟的时间步数。结合实验数据与COMSOL模型，进行联合仿真编程，代码如下：

functionResT = ObjFunc(x)
global xbest Fobjbest
tf          = 2500;
Cpel        = x(1)*1e6;
Cpcc        = x(1)*1e6;
keff        = x(2);
h1          = x(3);
ku        = x(4);
keref       = x(5);
alfa_ke     = x(6);
i0ref       = x(7);
Ei0         = x(8);
DS          = -0.07; 
Vinit       = 3.219;
%% Ambient Temperatures for Different Experiments
Tinf     = 23.8495;
%% Run COMSOL Model
LIT_A123_v1_0(Tinf,Cpel,Cpcc,keff,h1,ku,keref,alfa_ke,i0ref,Ei0,DS,Vinit)
%% Load EXPERIMENTAL Data
load SRun4c.mat
RunM(tf+1:end,:) = [];
tData            = RunM(:,1);
Tavg_data        = RunM(:,3);
Tmax_data        = RunM(:,2);
Tmin_data        = RunM(:,4);
Vcell_data       = RunM(:,5);
%% Load SIMULATION Results
Run1 = load('Run_4c_100s.txt');
close all
Run1(1,:) = [];
Run1(:,1) = Run1(:,1)-1; 
tModel      = Run1(:,1);
Tavg_mod    = Run1(:,3);
Tmax_mod    = Run1(:,2);
Tmin_mod    = Run1(:,4);
VcellModel  = Run1(:,5);
%% Calculate RESIDUALS
Res1        = (Tavg_mod-Tavg_data)./Tavg_mod;
Res2        = (Tmax_mod-Tmax_data)./Tmax_mod;
Res3        = (Tmin_mod-Tmin_data)./Tmin_mod;
Res4        = (VcellModel-Vcell_data)./VcellModel;
Res         = vertcat(Res1,Res2,Res3,Res4);
ResT        = Res(:);
Fobj        = sum(ResT.^2);
if Fobj<Fobjbest
    xbest    = x;
    Fobjbest = Fobj;
    save('last_res','xbest','Fobjbest')
end
%% PLOT EXP and SIM Results
figure('Position', [10, 1000, 600, 600])
plot(tData,Tavg_data,'-b',tData,Tmax_data,'-b',tData,Tmin_data,'-b',...
    tModel,Run1(:,2:end-2),'-g','linewidth',2)
ylabel('Temp(C)','fontsize',16)
xlim([0 tf])
title(['x =' num2str(x)],'fontsize',16)
figure('Position', [10, 0, 600, 600])
plot(tModel,VcellModel,'-g',tData,Vcell_data,'-b','linewidth',2)
xlabel('Time(sec)','fontsize',16)
ylabel('Voltage(V)','fontsize',16)
xlim([0 tf])
title(['Res = ' num2str(Fobj)],'fontsize',16)
x
Fobj
pause(0.1)