5
$\begingroup$

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

$\endgroup$
  • $\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
2
$\begingroup$

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.

$\endgroup$
  • $\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

protected by Community Aug 19 '15 at 17:16

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?

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