matlab - 使用有限差分求解 PDE

Question

我正在尝试使用有限差分法解决热失控反应的耦合 PDE。我有 2 个变量，温度 (T) 和浓度 (Cbulk)，它们随时间 (t) 和距离 (x) 的变化而变化。我已经为一个变量 (T) 解决了这个偏微分方程，同时保持 Cbulk 不变。

但是，当我尝试以温度和浓度作为变量来求解方程时，我得到了一个不合理的答案。在求解温度和浓度时，我附上了温度与 x 的关系图。谁能建议我能做些什么来纠正这个问题？实际答案应该是平滑的温度和浓度曲线。

另请参阅我的 MATLAB 代码

% fd_rhs.m 
function ydot=fd_rhs(time,y) 
global nt R delH E A rho cp h k delr Cbulk s

T = y(1:2:2*nt);
c = y(2:2:2*nt);


for i=1:nt-1
if i == 1
       dT(i) = k/(rho*cp*4.184*1000)*2*(s+1)*(T(i+1)-   
       T(i))/(delr*delr)+delH*A*exp(-E*1000./(R*T(i)))*c(i)/(rho*cp*1000);
       dc(i)= -A*exp(-E*1000./(R*T(i)))*c(i);
else 
dT(i) = k/(rho*cp*4.184*1000)*((T(i+1)-2*T(i) +T(i-1))/(delr*delr))+(s/R)* 
        (T(i+1)-T(i-1))/(delr*delr)+delH*A*exp(- 
        E*1000./(R*T(i)))*Cbulk/(rho*cp*1000); 
dc(i)= -A*exp(-E*1000./(R*T(i)))*c(i); 
end 
end

dT(nt)= (k/delr)/(h+(k/delr))*(k/(rho*cp*4.184*1000)*((T(nt)-2*T(nt-1) 
        +T(nt-2))/(delr*delr)))+(s/R)*(T(nt)-T(nt-1))/(delr*delr);
dc(nt)= -A*exp(-E*1000./(R*T(i)))*c(nt); 
ydot = [dT';dc']
end

谢谢

标签： matlabpde