# How to add reaction and source terms to a diffusion PDE solver written with MATLAB's pdepe function?

I have the following system of equations which I'm trying to solve using Matlab's pdepe solver.

The 1-D spherical heat diffusion equation with heat generation (source term):

$$\rho \, C_p\frac{\partial T}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r} \left(r^2 \, k \frac{\partial T}{\partial r} \right) + \Delta H \frac{\partial \rho_w}{\partial t} \\$$

and boundary conditions:

$$k \frac{\partial T}{\partial r} = h\,(T_\infty - T_s) \; \; \text{at r = surface}\\ \frac{\partial T}{\partial r} = 0 \; \; \text{at r = 0}$$

also note that $\Delta H$ is a constant, i.e. $\Delta H = -255,000$. The initial value for $\rho_w = \rho_{initial} = 700$ and the ambient temperature $T_\infty = 773$ (see Matlab code below for more details).

The kinetic reactions of the system:

$$\frac{\partial \rho_w}{\partial t} = -(K_1+K_2+K_3) \rho_w \\ \frac{\partial \rho_g}{\partial t} = K_1 \rho_w \\ \frac{\partial \rho_T}{\partial t} = K_2 \rho_w \\ \frac{\partial \rho_c}{\partial t} = K_3 \rho_w \\ \frac{\partial \rho_g}{\partial t} = K_4 \rho_T \\ \frac{\partial \rho_c}{\partial t} = K_5 \rho_T \\ \text{where}\; \; K = A\,e^{\frac{-E}{RT}}$$

The thermal parameters $\rho, C_p, k$ vary according to the following:

$$\rho = \rho_w + \rho_c \\ Y_w = \frac{\rho_w}{\rho_{initial}} \\ C_p = Y_w C_{pw} + (1-Y_w) C_{pc} \\ k = Y_w k_w + (1-Y_w) k_c$$

where $C_{pw} = 1500$, $C_{pc} = 1100$, $k_w = 0.105$, and $k_c = 0.071$.

I can solve the single PDE for heat diffusion (no source term) using the pdepe function in Matlab as follows:

function PDEexampleSphere

d350 = 0.035e-2; % diameter
Ti = 300;
Tinf = 773;
tmax = 0.8;

m = 2;
x = linspace(0,r350,20);
t = 0:0.01:tmax;

sol = pdepe(m,@pdefunc,@icfunc,@bcfunc,x,t);
u = sol;

% surface plot
figure(3)
surf(x,t,u)
xlabel('Distance (m)')
ylabel('Time (s)')

% temperature profile
figure(4)
plot(t,u(:,1),'b',t,u(:,end),'r')
hold on
plot([0 tmax],[Tinf Tinf],':k')
hold off
axis([0 tmax Ti-20 Tinf+20])
xlabel('Time (s)')
ylabel('Temperature (K)')

% --------------------------------------------------------------
function [c,f,s] = pdefunc(x,t,u,dudx)
rho = 700;  % density
cp = 1500;  % heat capacity
k = 0.105;  % thermal conductivity

c = rho*cp;
f = k.*dudx;
s = 0;

% --------------------------------------------------------------
function u0 = icfunc(x)
Ti = 300;   % initial temperature

u0 = Ti;

% --------------------------------------------------------------
function [pl,ql,pr,qr] = bcfunc(xl,ul,xr,ur,t)
h = 375;    % heat transfer coefficient
Tinf = 773; % ambient temperature

pl = 0;
ql = 0;
pr = -h*(Tinf-ur);
qr = 1;


How can I include the kinetic reactions and heat generation term into the Matlab function?

Note - for more information on the pdepe function in Matlab, please see the following documentation: http://www.mathworks.com/help/matlab/ref/pdepe.html

• This might be relevant... youtube.com/watch?v=ri_1nxwupb8 Oct 13, 2013 at 13:48
• @boyfarrell The video demonstrates exactly what I did in the above Matlab code for a single PDE and no source term. While helpful, it does not answer my question. Oct 13, 2013 at 14:26
• What do the docs say about adding a source term? Can't you literally just add a source term to your equation, is there something preventing you from doing that? Oct 13, 2013 at 23:39
• @boyfarrell The documentation states $s(x,t,u,\frac{\partial u}{\partial x})$ but my source term is in terms of $\rho$ not temperature. Oct 14, 2013 at 2:12
• Does the pdepe solver accept coupled PDEs? If so you could solve the first equation for $\frac{\partial T}{\partial t}$ and the second for $\frac{\partial u}{\partial t}$? You could simply use the solution of the second equation as the source term to the first. Provided you have initial conditions (so you know $s(t=0)$) it should be possible to include the value numerically in this way. Sorry I'm not really at MATLAB person, these are the things I would try. Oct 14, 2013 at 7:20

Since pdepe accepts systems of PDEs through vector-valued capacity, flux, and source terms, one way to accommodate your request would be to set the fluxes for all of the $\rho$ variables equal to zero. The capacity terms for the $\rho$ variables will all be 1, and the source terms for each variable are the non-flux terms on the right-hand side (the $\rho$ variables have no flux terms, so for these variables, the source term is also the only term on the right-hand side of these equations). Without knowing more about the methods used by pdepe, it's difficult to say whether or not this kludge will be numerically stable, but my guess is that it would probably work.