MATLAB锂离子电池P2D电化学仿真：从PDE离散到参数校准

简介：一套基于MATLAB的伪二维（P2D）模型仿真实现，即多伊尔-富勒-纽曼（DFN）模型，面向锂电池研发、电化学仿真方向的研究人员与工程师。其将电池复杂结构简化为一维电极厚度方向并引入粒子径向的伪第二维度，可同时模拟浓度与电势的时空演变，洞察比等效电路模型和单粒子模型更丰富的内部动力学。压缩包共130个文件，以27个.m脚本为主，覆盖电池参数定义、几何建模、有限元矩阵组装、P2D求解器配置与结果处理；另含100个.xml配置文件，以及.md说明文档、.mlx实时脚本和.prj工程文件，便于直接运行与二次开发。包体仅430KB，轻量紧凑，已有157人学习浏览。借助该代码，可快速复现P2D模型求解流程，调整材料参数和工况开展对比实验，并可视化电压、浓度等结果，为电池性能预测与优化提供实用平台。

1. 为什么要在 MATLAB 里做锂离子电池的 P2D 电化学仿真

做新能源项目时，最常遇到的建模分歧不是“用不用机器学习”，而是要不要从等效电路换成电化学模型。等效电路能快速拟合端电压，却回答不了“负极表面嵌锂量有没有越界”和“倍率上升时液相极化占大头还是固相极化占大头”，这两个问题都要靠 P2D 模型（pseudo-two-dimensional，准二维模型）在 MATLAB 里把锂离子电池的电化学行为重算一遍。P2D 沿电池厚度方向切一维，再把每颗活性颗粒沿半径方向切一维，用偏微分方程组描述固相扩散、液相扩散迁移和 Butler-Volmer 反应动力学。适合电芯设计、BMS 算法、储能系统工程师，用于复现充放电曲线、估算内阻和确认安全边界。下面按“机理—离散化—求解器—后处理”推进，用脚本路径跑通一个可解释的 P2D 模型。

2. P2D 模型的电化学机理与 MATLAB 离散化选型

2.1 从 PDE 方程组看 P2D 的核心变量与边界条件

P2D 并不是严格的二维模型，而是两个一维空间坐标耦合：电池厚度方向 (x) 和电极颗粒半径方向 (r)。在厚度方向上，依次是负极集流体、负电极、隔膜、正电极、正极集流体；在颗粒半径方向上，锂离子从颗粒表面向中心扩散。一个完整的 P2D 描述这些变量：固相锂浓度 (c_s(x,r,t))、液相锂浓度 (c_e(x,t))、固相电势 (\phi_s(x,t))、液相电势 (\phi_e(x,t))、局部反应电流密度 (j(x,t))。

控制方程可以简写成四组：固相菲克扩散方程、液相物料守恒方程、固相电荷守恒方程、液相电荷守恒方程。流量在界面处用 Butler-Volmer 方程耦合。求解边界条件不是随便取的，负极和正极集流体边界的固相电位梯度等于外加电流，隔膜两侧则要求液相通量连续，颗粒中心处固相浓度梯度为零。初学者最容易漏的，正是负极集流体边界和隔膜/电极界面的通量匹配。

2.2 为什么用有限差分而不是 MATLAB 的 PDE 工具箱

MATLAB 自带pdepe和 Partial Differential Equation Toolbox，可以解常见一维或二维抛物方程，但 P2D 是耦合多种物理场的微分代数方程：液相电势没有显式的时间导数项，pdepe并不擅长；另一个问题是固相颗粒内部半径方向和电池厚度方向需要同时离散，而两个方向的特征长度相差好几个数量级。常见做法是用有限差分做空间离散，把偏微分方程转换成常微分方程组，再用ode15s在时间维度上推进。这种“方法线”思路让每个物理量对应一个向量，调试和扩展副反应模型都更直接。

网格生成可以用 MATLAB 的向量化写法，下面这段代码把负极、隔膜、正极按实际厚度拼接成一套均匀网格：

% p2d_grid.m : 一维厚度方向网格生成 Nn = 20; Ns = 10; Np = 20; % 各区域网格数 Ln = 92e-6; Ls = 25e-6; Lp = 75e-6; % 负极/隔膜/正极厚度，单位 m % 先做全长度上的等距网格，再标记区域 Nx = Nn + Ns + Np; x_boundary = linspace(0, Ln + Ls + Lp, Nx + 1)'; x_mid = 0.5 * (x_boundary(1:end-1) + x_boundary(2:end)); % 用区域编号标记每个体中心：1=负极，2=隔膜，3=正极 region = zeros(Nx, 1); region(x_mid <= Ln) = 1; region(x_mid > Ln & x_mid <= Ln + Ls) = 2; region(x_mid > Ln + Ls) = 3;

参数说明：Nn/Ns/Np是三个区域的网格份数，不是节点个数；网格越密，浓度梯度和界面通量的解析越准，但状态数会同步增加。x_boundary存的是每个控制体的交界面，x_mid存的是控制体中心，实际离散液相浓度和电势都定义在体中心。region向量在做系数矩阵装配时非常关键，它可以避免后面反复用if-else判断当前网格属于哪个区域。

如果电极/隔膜界面上浓度梯度很大，等距网格会造成振荡，更稳的做法是在界面附近加密。用一个余弦函数映射就可以满足多数仿真需求：先在线性坐标zeta = linspace(0, pi, M+1)上等分，再取x = L * (1 - cos(zeta)) / 2，这样两端网格密、中间网格松。需要注意，加密后的网格要重新计算体心和面通量，不能用等距网格的中心差分系数。

2.3 空间尺度、时间尺度与收敛前的三条检查

P2D 模型最容易被忽略的是量纲。锂离子电池中浓度量级在 (10^3) 到 (10^4\ \mathrm{mol/m^3})，电势量级在 (0) 到 (5\ \mathrm{V})，而扩散系数量级在 (10^{-14}) 到 (10^{-10})，直接把这些数塞进ode15s，绝对误差容限很难设置。常见做法是在进入求解器前先归一化：把浓度除以初始电解液浓度或最大固相浓度，把空间坐标除以相应区域厚度，把时间除以特征扩散时间。这样所有状态量都在 1 附近，AbsTol可以设置成1e-6而不会把高量级大数误判为未收敛。

时间尺度也需要预判。固相扩散的时间常数通常是几百秒到上千秒，而液相浓度和电势几乎瞬时调整；整组方程是刚性系统，这也是使用ode15s而不是ode45的原因。如果ode15s每步都要做很多次迭代，先检查是否有代数变量没有正确表示成质量矩阵零行，再看网格长宽比是否超过 (10^6)。还有一个常见问题：初始固相浓度和液相浓度不一致时，第一时刻的 Butler-Volmer 过电势会被拉得很大，导致求解器在开始几微秒内疯狂缩短步长。解决办法是让初始浓度分布满足电荷中性条件，再缓慢加载电流，而不是直接给一个大电流阶跃。

提示：P2D 离散化完成后，先跑一个零电流静置工况。若零电流下端电压还会漂移，说明初始浓度或平衡电极电位不一致，先修物理参数，不要急着调求解器。

3. 在 MATLAB 中搭建 P2D 求解器的分步实现与参数

3.1 参数脚本：把物性参数整理成结构体

P2D 的参数量很大，不建议散落在私有函数里，用 MATLAB 结构体统一管理是最常见的做法。下面是一个最小参数集，后续所有函数都从p中读取：

% p2d_params.m : 返回参数结构体 function p = p2d_params() p.F = 96485; % 法拉第常数，C/mol p.Rg = 8.314; % 气体常数，J/(mol*K) p.T = 298.15; % 温度，K % 几何 p.Ln = 92e-6; p.Ls = 25e-6; p.Lp = 75e-6; p.Rsn = 6e-6; p.Rsp = 4e-6; % 颗粒半径，m % 液相 p.ce0 = 1000; % 初始电解液浓度，mol/m^3 p.De_ref = 1.0e-10; % 液相扩散参考值，m^2/s p.brug = 1.5; % Bruggeman 曲折因子指数 p.tplus = 0.363; % 阳离子迁移数 % 固相 p.csmax_n = 31507; % 负极最大嵌锂浓度，mol/m^3 p.csmax_p = 22806; % 正极最大嵌锂浓度，mol/m^3 p.Ds_n = 3.0e-14; % 负极固相扩散系数，m^2/s p.Ds_p = 3.7e-14; % 正极固相扩散系数，m^2/s % 反应动力学参考交换电流密度，A/m^2 p.k0_n = 2.0e-11; p.k0_p = 3.0e-11; end

参数说明：液相扩散系数和交换电流密度是温度相关的，这里的De_ref和k0_n/k0_p只是常温参考值，实际工程中应该用 Arrhenius 公式修正温度。brug通常取 1.5，代表多孔电极的曲折效应；提高brug等效于降低液相有效导电率，会让倍率放电时的液相极化明显变大。tplus一般小于 0.5，它直接进入液相物料守恒方程中的对流迁移项，这个值错 0.1 就可能让高倍率下的浓度曲线整体偏离。

3.2 组装右端项函数：固相扩散与液相扩散的离散骨架

求解 P2D 的代码结构通常是一个主脚本加一个残差函数。残差函数的输入是归一化状态向量x，输出是dxdt；电势和反应电流密度在每个时间步内通过代数方程解出。为了把代码控制在可维护范围内，固相径向和液相厚度方向的二阶导都预先生成稀疏差分矩阵：

% build_matrices.m : 生成固相径向、液相厚度方向的二阶差分矩阵 function [A_s, A_e] = build_matrices(p, mesh) Nr = mesh.Nr; Nx = mesh.Nx; % 固相颗粒半径方向：中心差分，r = [0, Rs] dr = p.Rsn / (Nr - 1); r = (0:Nr-1)' * dr; A_s = spdiags([...], -1:1, Nr, Nr); % 实际需按边界修正 r=0 和 r=Rs 项 % 液相厚度方向：有效扩散系数随区域变化，先由 region 生成系数向量 Deff = p.De_ref * (mesh.porosity .^ p.brug); % A_e 使用通量守恒型离散，保证电极/隔膜界面上的连续性 A_e = assemble_conservative_laplacian(Deff, mesh.dx); end

逻辑说明：固相径向采用球坐标菲克扩散，所以一阶导数的系数会带 (1/r^2)，如果没有在r=0处做边界处理，矩阵会奇异。液相厚度方向优先使用通量守恒型离散，也就是先计算界面通量，再减控制体两侧通量差；如果直接用标准五点中心差分，隔膜和电极界面的有效扩散系数跳变会产生虚假通量。

实际操作中不必每次重新生成差分矩阵，可以把A_s和A_e做成稀疏矩阵，在ode15s调用之前只生成一次。注意 MATLAB 的spdiags在处理变系数时不宜直接传入向量，推荐用diag或自定义装配函数把离散项叠加到稀疏矩阵上。代码里的assemble_conservative_laplacian就是用来处理变系数拉普拉斯算子的自定义函数，它返回一个 (N_x \times N_x) 的稀疏矩阵，运行时间远小于读入参数的时间。

3.3 用 ode15s 求解刚性 P2D 方程组的调用命令

把状态向量组织成[固相浓度向量; 液相浓度向量]后，主程序调用ode15s。这里需要传入质量矩阵，把代数方程对应的行设成 0：

% run_p2d.m : 主求解入口 [t, x] = ode15s(@(t, x) p2d_rhs(t, x, p, mesh), ... [0, 3600], x0, options); % 输出端电压：每个时刻都要先解出表面过电势和膜电阻压降 V = zeros(size(t)); for k = 1:numel(t) V(k) = compute_voltage(t(k), x(k,:).', p, mesh); end

options通常是：

options = odeset('RelTol', 1e-4, 'AbsTol', 1e-6, ... 'Mass', M, 'MassSingular', 'yes', ... 'Jacobian', @(t, x) p2d_jac(t, x, p, mesh));

参数说明：M是质量矩阵，微分变量所在行对角线为 1，代数变量所在行为 0；MassSingular设置为'yes'，ode15s会按索引 1 的微分代数方程处理。Jacobian可选但推荐提供，P2D 的状态数在 500 到 2000 之间，不提供解析雅可比时ode15s会做数值差分，每次迭代额外开销很大，而且容易因状态量量级差异产生截断误差。若不想手推雅可比，至少把JPattern设成稀疏结构矩阵，让ode15s只对非零元素做数值差分。

倍率变化不要直接写成一个大幅值阶跃电流，常见做法是把电流序列I(t)线性插值进残差函数。这样每个步长都能得到连续电流，ode15s不会在电流跳变点反复缩短步长。compute_voltage要包含热力学平衡电位、固相和液相的过电势以及 SEI 膜电阻压降，不能只看 Butler-Volmer 算出来的反应过电势，否则端电压曲线会和实测差出几十毫伏。

4. 用 P2D 模型分析电化学行为：倍率、内阻与敏感性

4.1 倍率扫描：从端电压曲线看极化分配

P2D 模型最有价值的输出不是“电压多准确”，而是能把端电压拆开看。在 MATLAB 中做倍率扫描非常简单，对 1C、3C、5C 分别调用一次求解器，然后画在同一张图上：

C_rate = [1, 3, 5]; colors = lines(3); figure('Color', 'w'); hold on; for k = 1:numel(C_rate) I_app = C_rate(k) * p.capacity; % p.capacity 为电芯容量 A*h [t, x] = run_p2d(p, mesh, I_app); V = compute_voltage(t, x, p, mesh); plot(t / 3600, V, 'LineWidth', 1.6, 'Color', colors(k, :)); end xlabel('时间 (h)'); ylabel('端电压 (V)'); legend('1C', '3C', '5C'); grid on;

逻辑说明：p.capacity由电极活性物质体积、最大浓度和初始嵌锂程度共同决定，不能只填任意整数；正确做法是从参数化模型中计算出理论容量，再乘以库仑效率系数。倍率越大，早期压降越陡，原因是欧姆极化瞬间建立，而浓差极化需要一段时间才充分发展。把 1C 和 5C 的电压差画成另一条曲线，就可以看到“液相扩散慢”在后期贡献的增长趋势，这是用等效电路很难直接分离出来的信息。

对于 5 年以上经验的工程师，更建议在倍率扫描后同时输出负极表面固相浓度 (c_{s,\mathrm{surf}}) 的时间曲线。P2D 里真正决定安全边界的不是端电压，而是负极表面嵌锂分数。如果表面浓度在大倍率放电末端逼近最大浓度，即使端电压还没到截止电压，也应该在 BMS 策略里提前限功率。

4.2 参数敏感性：先扰动后扫描的 MATLAB 操作

研究 P2D 模型的电化学行为，不能只调一个参数看一条曲线。推荐先做单参数扰动，找到响应梯度大的参数，再做小范围扫描。下面这段代码用一个两层循环做批量求解，把最终放电容量记录成表格：

% sensitivity_sweep.m param_list = {'Ds_n', 'De_ref', 'brug', 'k0_p'}; ratio = [0.5, 0.75, 1.0, 1.25, 1.5]; S = table(); for pi = 1:numel(param_list) for ri = 1:numel(ratio) p_tmp = p; p_tmp.(param_list{pi}) = p.(param_list{pi}) * ratio(ri); cap = discharge_capacity(p_tmp, mesh, I_app); S = [S; table(param_list{pi}, ratio(ri), cap)]; end end

参数说明：discharge_capacity需要自己封装求解过程和电压截止判断，返回值为放电至截止电压时释放的容量。S用表类型存储，方便后续直接画groupedbar或写入 Excel。装配参数矩阵时不要用eval，MATLAB 结构体字段动态访问用p_tmp.(param_list{pi})既快又安全。

从常见文献和工程数值看，倍率工况下De_ref和brug的敏感性通常高于k0，因为高倍率时液相传输成为主要瓶颈；低倍率工况下k0和csmax的影响更明显。这个结论不是绝对的，颗粒半径Rsn/Rsp改变也会移动瓶颈位置。如果批量仿真时间太长，可以先用 Simulink 或并行parfor把参数扫描拆到多核上；注意parfor里不要动态修改同一个结构体p的字段，要先拷贝一份p_tmp。

表格里的“响应”指端电压或容量的变化幅度。实际项目里修改任何一个参数，都必须同步检查电池的 OCV 曲线是否改变；很多参数只应该影响动力学项，不应改变平衡电位，如果模型把平衡电位写成定值而不是浓度插值函数，敏感性分析的结果会失真。

4.3 内阻贡献拆解：从求解结果重新合成端电压

端电压在放电过程中会同时包含欧姆压降、活化极化、浓差极化。P2D 的好处是每个贡献项都对应明确的物理项。可以在残差函数里把每个时刻的液相过电势和固相过电势分别缓存到求解器输出结构中，然后拆开画：

% 在 compute_voltage 中额外返回分量 [V, eta_ohm, eta_act, eta_conc] = compute_voltage(t, x, p, mesh); figure; area(t / 3600, [eta_ohm, eta_act, eta_conc] .* I_app); legend('欧姆', '活化', '浓差');

eta_ohm来自液相离子导电率和固相电子导电率产生的电位梯度，eta_act来自 Butler-Volmer 方程的交换电流密度，eta_conc来自颗粒表面和体相浓度的浓度差。这样拆开后，如果某款电芯高倍率下浓差极化占比超过 50%，多半要优先改电解液电导率或电极厚度，而不是单纯调 SEI 膜电阻。这个“拆积木”的能力是电池等效电路模型很难做到的，也是 P2D 在 MATLAB 中最实用的电化学行为分析手段。

5. 把 P2D 仿真结果变成可验证数据的三个技巧

5.1 用插值表校准平衡电极电位，不要用常数

很多 P2D 入门代码把正负极平衡电位设成固定值，导致放电平台和实测台阶对不上。正确做法是把正负极的嵌锂量转为 SOC，再用 MATLABinterp1查表：

Eeq_n = interp1(soc_n_table, ocp_n_table, soc_n, 'pchip'); Eeq_p = interp1(soc_p_table, ocp_p_table, soc_p, 'pchip');

这里的soc_n用当前固相表面浓度除以最大浓度得到；ocp_n_table要从半电池实验获得，不能用成品电芯实测曲线替代。拟合后要检查正负极电位曲线是否在循环过程中越过平台区，如果越过了，说明模型容量或初始嵌锂量设定错了。

5.2 用多倍率数据按“先固相后液相”的顺序校准

校准参数时有固定的先后顺序：先用低倍率找平衡电位和初始嵌锂量，再用中倍率校准固相扩散系数Ds_n/Ds_p，最后用高倍率校准液相参数De_ref/brug。原因很简单，低倍率下液相浓差还没有建立，端电压主要由热力学和活化极化控制；高倍率下液相扩散表现最突出。这个顺序反了，容易得到两个参数互相补偿的错误组合。

5.3 用事件函数在负极析锂边界自动切断放电

ode15s支持事件函数，可以在负极表面嵌锂浓度到达阈值时自动终止仿真。这个技巧比写死截止电压更接近真实保护逻辑：

function [value, isterminal, direction] = pli_event(~, x, p, mesh) cs_surf = extract_negative_surface_concentration(x, p, mesh); value = p.csmax_n * 0.98 - cs_surf; % 降到 0 时触发 isterminal = 1; direction = 0; end

value是阈值条件，降到 0 表示负极表面浓度已到上限；isterminal = 1让求解器停止；direction = 0表示正负方向穿过零都触发。设置0.98是给模型误差留余量，实际值由电芯厂商给出。这样跑完一次放电，t(end)和最后的电压就是“该工况下能够安全放出的电量”，也是 BMS 策略里非常有价值的仿真输入。

把这三件事做成通用的前处理和后处理函数，P2D 模型的仿真结果就能和实验室倍率曲线、半电池 OCV 数据直接对齐，后续加副反应或析锂模型时，也不需要推翻原有求解框架。

本文还有配套的精品资源，点击获取