🎯要点
- 模拟对比石墨电池的放电电压曲线与实验数据定性差异。
- 对比双箔、多相多孔电极理论和锂电有限体积模型实现。
- 通过孔隙电极理论模型了解粗粒平均质量和电荷传输以及孔隙率的表征意义。
- 锂电中锂离子正向和逆向反应速率与驱动力的指数以及电解质和电极表面的锂浓度在经验上相关。
- 固相物质传递和电解质相中的质量传输过程的偏微分方程。
🍁电池电化学和热力学
🍪语言内容分比
🍇MATLAB偏微分方程
一维偏微分方程包含一个依赖于时间 t t t 和一个空间变量 x x x 的函数 u ( x , t ) u(x, t) u(x,t)。MATLAB PDE 求解器 pdepe 可求解以下形式的一维抛物线和椭圆 PDE 系统,
c ( x , t , u , ∂ u ∂ x ) ∂ u ∂ t = x − m ∂ ∂ x ( x m f ( x , t , u , ∂ u ∂ x ) ) + s ( x , t , u , ∂ u ∂ x ) c\left(x, t, u, \frac{\partial u}{\partial x}\right) \frac{\partial u}{\partial t}=x^{-m} \frac{\partial}{\partial x}\left(x^m f\left(x, t, u, \frac{\partial u}{\partial x}\right)\right)+s\left(x, t, u, \frac{\partial u}{\partial x}\right) c(x,t,u,∂x∂u)∂t∂u=x−m∂x∂(xmf(x,t,u,∂x∂u))+s(x,t,u,∂x∂u)
偏导数相对于时间的耦合仅限于乘以对角矩阵 c ( x , t , u , ∂ u ∂ x ) c\left(x, t, u, \frac{\partial u}{\partial x}\right) c(x,t,u,∂x∂u)。此矩阵的对角线元素要么为零,要么为正数。零元素对应于椭圆方程,其他元素对应于抛物线方程。必须至少有一个抛物线方程。如果 c c c 中对应于抛物线方程的元素是网格点(求解的点),则它们可以在孤立的 x x x 值处消失。只要在每个界面上放置一个网格点,就可以允许因材料界面而导致 c c c 和 s s s 不连续。
要使用 pdepe 求解 PDE,您必须定义 c 、 f c、f c、f 和 s s s 的方程系数、初始条件、解在边界上的行为以及用于计算解的点网格。函数调用 sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)
使用此信息在指定的网格上计算解:
必须以 pdepe 所要求的标准形式来表达偏微分方程。以这种形式书写,可以读出系数 c 、 f c、f c、f 和 s s s 的值。在 MATLAB 中,可以用以下形式的函数对方程进行编码:
function [c,f,s] = pdefun(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end
在这种情况下,pdefun 定义方程 ∂ u ∂ t = ∂ 2 u ∂ x 2 \frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2} ∂t∂u=∂x2∂2u。如果有多个方程,则 c 、 f c、f c、f 和 s s s 是向量,每个元素对应一个方程。在初始时间 t = t 0 t=t_0 t=t0,对于所有 x x x,解的分量满足如下形式的初始条件 u ( x , t 0 ) = u 0 ( x ) u\left(x, t_0\right)=u_0(x) u(x,t0)=u0(x).
在 MATLAB 中,可以用以下形式的函数对初始条件进行编码
function u0 = icfun(x)
u0 = 1;
end
在这种情况下, u ∅ = 1 u \varnothing=1 u∅=1 定义了 u 0 ( x , t 0 ) = 1 u_0\left(x, t_0\right)=1 u0(x,t0)=1 的初始条件。如果有多个方程,则 u θ u \theta uθ 是一个向量,每个元素定义一个方程的初始条件。
在边界 x = a x=a x=a 或 x = b x=b x=b 处,对于所有 t t t,解的各分量均满足以下形式的边界条件
p ( x , t , u ) + q ( x , t ) f ( x , t , u , ∂ u ∂ x ) = 0 p(x, t, u)+q(x, t) f\left(x, t, u, \frac{\partial u}{\partial x}\right)=0 p(x,t,u)+q(x,t)f(x,t,u,∂x∂u)=0
q ( x , t ) q(x, t) q(x,t) 是一个对角矩阵,其元素要么为零,要么不为零。请注意,边界条件以通量 f f f 表示,而不是 u u u 相对于 x x x 的偏导数。此外,在两个系数 p ( x , t , u ) p(x, t, u) p(x,t,u) 和 q ( x , t ) q(x, t) q(x,t) 中,只有 p p p 可以依赖于 u u u。
在 MATLAB 中,可以用以下形式的函数对边界条件进行编码
function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t)
pL = uL;
qL = 0;
pR = uR - 1;
qR = 0;
end
pL 和 qL 是左边界的系数,而 pR 和 qR 是右边界的系数。在这种情况下,bcfun 定义边界条件
u L ( x L , t ) = 0 u R ( x R , t ) = 1 \begin{aligned} & u_L\left(x_L, t\right)=0 \\ & u_R\left(x_R, t\right)=1 \end{aligned} uL(xL,t)=0uR(xR,t)=1
如果有多个方程,则输出 pL、qL、pR 和 qR 是向量,每个元素定义一个方程的边界条件。
MATLAB PDE 求解器中的默认积分属性用于处理常见问题。在某些情况下,您可以通过覆盖这些默认值来提高求解器性能。为此,使用 odeset 创建一个选项结构。然后,将该结构作为最后一个输入参数传递给 pdepe:sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)
抛物线型偏微分方程的一个例子是一维热方程:
∂ u ∂ t = ∂ 2 u ∂ x 2 \frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2} ∂t∂u=∂x2∂2u
在对方程进行编码之前,需要确保它符合 pdepe 求解器所期望的形式:
c ( x , t , u , ∂ u ∂ x ) ∂ u ∂ t = x − m ∂ ∂ x ( x m f ( x , t , u , ∂ u ∂ x ) ) + s ( x , t , u , ∂ u ∂ x ) c\left(x, t, u, \frac{\partial u}{\partial x}\right) \frac{\partial u}{\partial t}=x^{-m} \frac{\partial}{\partial x}\left(x^m f\left(x, t, u, \frac{\partial u}{\partial x}\right)\right)+s\left(x, t, u, \frac{\partial u}{\partial x}\right) c(x,t,u,∂x∂u)∂t∂u=x−m∂x∂(xmf(x,t,u,∂x∂u))+s(x,t,u,∂x∂u)
在此形式中,热方程为
1 ⋅ ∂ u ∂ t = x 0 ∂ ∂ x ( x 0 ∂ u ∂ x ) + 0 1 \cdot \frac{\partial u}{\partial t}=x^0 \frac{\partial}{\partial x}\left(x^0 \frac{\partial u}{\partial x}\right)+0 1⋅∂t∂u=x0∂x∂(x0∂x∂u)+0
因此系数的值如下:
- m = 0 m=0 m=0
- c = 1 c=1 c=1
- f = ∂ u ∂ x f=\frac{\partial u}{\partial x} f=∂x∂u
- s = 0 s=0 s=0
使用 20 个点的空间网格和 30 个点的时间网格。由于解迅速达到稳定状态,因此 t = 0 t=0 t=0 附近的时间点间隔更紧密,以便在输出中捕捉这种行为。
L = 1;
x = linspace(0,L,20);
t = [linspace(0,0.05,20), linspace(0.5,5,10)];
最后,利用对称性 m m m、PDE方程、初始条件、边界条件以及 x x x 和 t t t 的网格求解该方程。
m = 0;
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);
使用 pcolor
可视化解矩阵。
colormap hot
pcolor(x,t,sol)
colorbar
xlabel('Distance x','interpreter','latex')
ylabel('Time t','interpreter','latex')
title('Heat Equation for $0 \le x \le 1$ and $0 \le t \le 5$','interpreter','latex')
局部函数
function [c,f,s] = heatpde(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end
function u0 = heatic(x)
u0 = 0.5;
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur - 1;
qr = 0;
end