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.

enter image description here

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)-   
   dc(i)= -A*exp(-E*1000./(R*T(i)))*c(i);
dT(i) = k/(rho*cp*4.184*1000)*((T(i+1)-2*T(i) +T(i-1))/(delr*delr))+(s/R)* 
dc(i)= -A*exp(-E*1000./(R*T(i)))*c(i); 

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

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.

Temperature vs distance

Concentration vs distance

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.Temperature vs distance at very small time = 0.0001 hours

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

  • $\begingroup$ 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? $\endgroup$ – Anton Menshov Dec 27 '18 at 16:39
  • $\begingroup$ 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. $\endgroup$ – user10837299 Dec 27 '18 at 16:45
  • 1
    $\begingroup$ 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? $\endgroup$ – Savithru Dec 27 '18 at 18:06
  • 1
    $\begingroup$ 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. $\endgroup$ – Bill Greene Dec 27 '18 at 20:05
  • $\begingroup$ Thank you for your guidance. As I read, I think PDEPE might be worth a try. Thank you. $\endgroup$ – user10837299 Dec 28 '18 at 2:39

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.