# Finite difference coupled PDE solution to 1D heat-reaction equation

I am trying to solve a coupled PDE for a thermal runaway reaction using finite difference method. I have 2 variables, temperature (T) and concentration (c) that vary as a function of time (t) and distance (x). I have already solved this PDE for one variable (T) while keeping c constant.

Now I am trying to solve the coupled PDE (both T and c).The equations are as follows.

My matlab code to solve this using ODE15s is as follows:

% 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)+H*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)+H*A*exp(-
E*1000./(R*T(i)))*c(i)/(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


When I try to solve the equation with both temperature and concentration as the variables, I am getting an irrational answer. I have attached a plot of temperature vs x when solving for both temp and concentration. Can anyone please suggest what I can do to correct this ? The actual answer should be a smooth temperature and concentration curve.

Even when I try to solve using small time steps (t=0.00001 hours) - there is a discontinuity that appears in the temperature curve at half the distance (at r=0.5 for R=1) - shown in the figure below It is always at half the distance -when R=2, it appears at r=1 and so on.

Thank you for your help and your time. I am hoping that someone could guide me in the right direction.

• Can you clarify the formula? You have both $C$ and $c$ there ($C_\text{bulk}$ also, but it is probably a constant). Are they different? or the same? What are $H$ and $A$ with a strange $.$ operator? – Anton Menshov Dec 27 '18 at 16:39
• Thank you for the comment. I have edited the code to make Cbulk as c(i). C bulk was the value of c when I held concentration constant to solve the PDE for only 1 variable(T). The H and A are the enthalpy of reaction and Pre-exp factor. The E is the reaction activation energy. These are all constant parameters in this PDE. – user10837299 Dec 27 '18 at 16:45
• Your MATLAB code doesn't seem to match the equations you have given. In your expressions for dT(i), you have a term that's multiplied by (s/R), but I don't see that term in the equations you have written. Similarly, you are missing the discretization of the $D\nabla^2 C$ term in your code for dc(i). Also, what are the boundary conditions for this problem? – Savithru Dec 27 '18 at 18:06
• Why don't you simply use the matlab pdepe function to solve these equations? It is very well-suited to solving systems of PDEs of this type that have time and a single spatial dimension as independent variables. – Bill Greene Dec 27 '18 at 20:05
• Thank you for your guidance. As I read, I think PDEPE might be worth a try. Thank you. – user10837299 Dec 28 '18 at 2:39