基于MATLAB求解Simulink金属氢化物吸放耦合PDE问题

一、可行性结论 #

完全可以在MATLAB环境中完成该耦合PDE系统的求解，pdepe工具是适配这类问题的核心工具，以下是具体适配思路与解决方案：

二、pdepe工具的适配方法 #

pdepe支持求解耦合抛物型/椭圆型PDE系统，你的问题需将3个耦合方程（2个质量平衡+1个能量平衡，耦合压力$P(r,t)$与温度$T(r,t)$）转换为其标准形式：

$$c(x,t,u,\frac{\partial u}{\partial x})\frac{\partial u}{\partial t} = x^{-m}\frac{\partial}{\partial x}\left(x^m f(x,t,u,\frac{\partial u}{\partial x})\right) + s(x,t,u,\frac{\partial u}{\partial x})$$
其中未知量向量$u = [P(r,t), T(r,t)]$，圆柱径向坐标系下$m=1$（对应权重项$r^1$）

1. 质量平衡方程的处理 #

无需改写含$\rho_g(P,T)$、$\rho_s(P,T)$的耦合项，只需在定义pdepe的核心函数时，将密度项作为$u(1)$（压力）和$u(2)$（温度）的函数直接代入。例如：

function [c,f,s] = mh_pde(r,t,u,dudr)
    % 提取未知量
    P = u(1);
    T = u(2);
    % 调用自定义函数计算密度等参数
    rho_g = calc_rhog(P,T); 
    rho_s = calc_rhos(P,T);
    % 整理第一个质量平衡方程的c、f、s项
    c(1) = ...; % 时间导数项系数
    f(1) = ...; % 扩散项（含dP/dr）
    s(1) = ...; % 源项（含P、T耦合项）
    % 第二个质量平衡方程同理
    c(2) = ...;
    f(2) = ...;
    s(2) = ...;
end

2. 能量平衡方程的适配 #

直接将能量方程拆解为pdepe要求的三部分：

• $c$：时间导数项的系数（如热容相关项）
• $f$：扩散项（如热导率与温度梯度的乘积）
• $s$：源项（如吸脱附反应热、对流换热项）

三、参考论文中方程转换的逻辑与求解 #

论文中将方程转换为无量纲形式（引入$\Theta$作为无量纲温度），核心目的是简化方程尺度、消除量纲干扰：

1. $\Theta$是无量纲温度，一般形式为$\Theta = \frac{T-T_{ref}}{T_{ref}}$，将实际温度转换为0-1区间的无量纲量
2. 转换后所有变量（压力、温度、径向坐标、时间）均为无量纲形式，方程系数更简洁，能提升数值求解稳定性
3. 求解时只需先定义无量纲参数，用pdepe求解得到无量纲的$\Theta$与压力后，再转换回有量纲量即可

四、具体实现步骤 #

1. 方程整理：将有量纲（或无量纲）方程严格对应pdepe的标准形式
2. 辅助函数编写：实现$\rho_g(P,T)$、$\rho_s(P,T)$、热导率、反应热等参数的计算函数
3. 初始/边界条件定义：
   • 初始条件函数：返回径向各点的初始压力与温度
   • 边界条件函数：对应中心对称条件、罐壁冷却换热条件等
4. 求解调用：设置径向网格r、时间点t，执行sol = pdepe(1, @mh_pde, @mh_ic, @mh_bc, r, t)（$m=1$对应圆柱坐标系）
5. 结果可视化：用mesh、surf等函数绘制$P(r,t)$与$T(r,t)$的时空分布

内容的提问来源于stack exchange，提问作者Jj no