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
xlabel('Distance (m)')
ylabel('Time (s)')

% temperature profile
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

  • $\begingroup$ This might be relevant... youtube.com/watch?v=ri_1nxwupb8 $\endgroup$
    – boyfarrell
    Oct 13 '13 at 13:48
  • $\begingroup$ @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. $\endgroup$
    – wigging
    Oct 13 '13 at 14:26
  • $\begingroup$ 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? $\endgroup$
    – boyfarrell
    Oct 13 '13 at 23:39
  • $\begingroup$ @boyfarrell The documentation states $s(x,t,u,\frac{\partial u}{\partial x})$ but my source term is in terms of $\rho$ not temperature. $\endgroup$
    – wigging
    Oct 14 '13 at 2:12
  • $\begingroup$ 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. $\endgroup$
    – boyfarrell
    Oct 14 '13 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.

  • $\begingroup$ Can you provide some Matlab code demonstrating this? $\endgroup$
    – wigging
    Oct 19 '13 at 3:47
  • $\begingroup$ Actually, some Python code demonstrating this would be better. $\endgroup$
    – wigging
    Aug 5 '16 at 14:24

Not the answer you're looking for? Browse other questions tagged or ask your own question.