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

Question

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
r350 = d350/2;   % radius
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

Answer 1

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.