I am using the implicit finite difference method to discretize the 1-D transient heat diffusion equation for solid spherical and cylindrical shapes:

$$ \frac{1}{\alpha}\frac{\partial T}{\partial t} = \frac{\partial^2 T}{\partial t^2} + \frac{p}{r} \frac{\partial T}{\partial r} \; \; \; \text{for} \; r\neq0 \\ \frac{1}{\alpha}\frac{\partial T}{\partial t} = (1+p)\frac{\partial^2 T}{\partial r^2} \; \; \; \text{for} \; r=0 \\ \text{note that }\; \; \alpha = \frac{k}{\rho C_p} $$

where $p=1$ for cylinder and $p=2$ for sphere.

The boundary conditions are: $$ \left.\begin{matrix} \frac{\partial T}{\partial r} \end{matrix}\right|_{r=0} = 0 \; \; \; \text{for center node}\\ \left.\begin{matrix} k\frac{\partial T}{\partial r} \end{matrix}\right|_{r=R} = h(T_\infty - T_S) \; \; \; \text{for surface node} $$ where $T_S$ is temperature at surface node and $R$ is outer radius of cylinder or sphere.

Using the above equations and boundary conditions I arrived at the following discretized approximations for the temperatures at radial points from the center to surface:

for the center node where $i=0$

$$ \left [ 1 + 2(1+p) Fo \right ]T_0^{\;n+1} - 2(1+p)FoT_1^{\;n+1} = T_0^\;n $$

for the internal nodes where $i=1,2,...,M-1$

$$ Fo\left ( 1-\frac{p}{2i} \right )T_{i-1}^{\;n+1}+(1+2Fo)T_i^{\;n+1}-Fo\left ( 1+\frac{p}{2i} \right )T_{i+1}^{\;n+1}=T_i^{\;n} $$

and finally for the surface node where $i=M$

$$ -2FoT_{M-1}^{\;n+1}+\left [ 1+2Fo \left (1+Bi+Bi\frac{p}{2M} \right ) \right ] T_M^{\;n+1} = T_M^{\;n} + 2FoBi\left ( 1+\frac{p}{2M} \right )T_\infty $$

where $n$ is the present time while $n+1$ in the next time level and $Fo=\alpha\Delta t / \Delta r^2$, $Bi=h\Delta r / k$, $\alpha = k/\rho c_p$, $h$ is heat transfer coefficient, $\rho$ is density, $k$ is thermal conductivity, $c_p$ is heat capacity.

So using the numerical equations, one can solve for the temperatures inside the sphere or cylinder by creating a system of equations in the form of $[A]\left \{ T \right \}=\left \{ B \right \}$ and solve the temperature at each node by using the Matlab operation T = A \ B

$$ \begin{bmatrix} 1+2(1+p)Fo & -2(1+p) & 0 & 0 & 0\\ Fo\left ( 1-\frac{p}{2i} \right ) & 1+2Fo & Fo\left ( 1+\frac{p}{2i} \right ) & 0 & 0\\ 0 & Fo\left ( 1-\frac{p}{2i} \right ) & 1+2Fo & Fo\left ( 1+\frac{p}{2i} \right ) & 0\\ 0 & 0 & Fo\left ( 1-\frac{p}{2i} \right ) & 1+2Fo & Fo\left ( 1+\frac{p}{2i} \right )\\ 0 & 0 & 0 & -2Fo & 1+2Fo \left (1+Bi+Bi\frac{p}{2M} \right ) \end{bmatrix} \begin{bmatrix} T_0^{\;n+1}\\ T_1^{\;n+1}\\ T_2^{\;n+1}\\ T_3^{\;n+1}\\ T_4^{\;n+1} \end{bmatrix} = \begin{bmatrix} T_0^{\;n}\\ T_1^{\;n}\\ T_2^{\;n}\\ T_3^{\;n}\\ T_4^{\;n}+2FoBi\left ( 1+\frac{p}{2M} \right )T_\infty \end{bmatrix} $$

But how do I include kinetic reactions into this system?

I have the following reaction rates for $w$ wood and $c$ char:

$$ \frac{\partial \rho_w}{\partial t} = -(K_1+K_2)\rho_w \\ \frac{\partial \rho_{c1}}{\partial t} = K_2\rho_w \\ \frac{\partial \rho_{c2}}{\partial t} = K_3\rho_{c1} $$

The rate constants $K$ are represented by the Arrhenius equation:

$$ K=A\,e^{\frac{-E}{RT}} $$

where $A$ = pre-factor, $E$ = activation energy, $R$ = universal gas constant, and $T$ = temperature at that node.

So to try to incorporate these reaction equations into my system of equations for the temperatures I have discretized the reactions using the implicit method:

$$ \left [ 1+(K_1+K_2)\Delta t \right ] \rho_{wi}^{\;n+1} = \rho_{wi}^{\;n} \\ \rho_{c1i}^{\;n+1}-K_2\rho_{wi}^{\;n+1}\Delta t = \rho_{c1i}^{\;n} \\ \rho_{c2i}^{\;n+1}-K_3\rho_{c1i}^{\;n+1}\Delta t = \rho_{c2i}^{\;n} $$

Any suggestion on how to incorporate these kinetic reactions into the system of temperature equations?

Or should I solve for the temperatures first, then use the new temperatures in the reaction rates to update the $\rho$, then use the updated $\rho$ for the next iteration?

  • $\begingroup$ The temperature equation you stated does not seem to depend on the $\rho$s. Given that, how do want these quantities to enter the temperature equation? $\endgroup$ Commented Sep 26, 2013 at 1:43
  • $\begingroup$ @WolfgangBangerth The heat diffusion equation depends on $\alpha$ which is the thermal diffusivity. As stated, $\alpha = \frac{k}{\rho c_p}$ so the thermal properties of $k$, $\rho$, and $c_p$ all effect the heat equation. $\endgroup$
    – wigging
    Commented Sep 26, 2013 at 2:35
  • $\begingroup$ @WolfgangBangerth Also note that in the numerical equations the term $Fo = \frac{\alpha \Delta t}{\Delta r^2}$ would deal with $\rho$, $k$, and $c_p$ $\endgroup$
    – wigging
    Commented Sep 26, 2013 at 3:24
  • $\begingroup$ And $\rho$ depends on $\rho_{c1},\rho_{c2},\rho_{wi}$? $\endgroup$ Commented Sep 27, 2013 at 2:00
  • $\begingroup$ @WolfgangBangerth Yes. The fraction of wood remaining at each time step at each node is $Yw = \rho_{wi} / \rho_w$ so the density at any time is actually $\rho = Yw*\rho_w + (1-Yw)*(\rho_{c1}+\rho_{c2})$. $\endgroup$
    – wigging
    Commented Sep 27, 2013 at 3:00

1 Answer 1


As I said in my answer to your previous question, it's probably better if you don't try to write your own ODE solver, which is what you are doing now. There are a lot of good libraries out there that will solve a system of ordinary differential equations much better than if you roll your own implementation. If you can use them, you should; since your previous post uses MATLAB, you have many different ODE solvers available to you, and there are additional ODE solvers you can install that improve upon these (such as SUNDIALS).

You'll have a system of equations that looks like:

\begin{align} \frac{\mathrm{d}T_{i}}{\mathrm{d}t} = \ldots, \\ \frac{\mathrm{d}w_{i}}{\mathrm{d}t} = \ldots, \\ \frac{\mathrm{d}c_{i}}{\mathrm{d}t} = \ldots, \\ \frac{\mathrm{d}a_{i}}{\mathrm{d}t} = \ldots, \\ \end{align}


  • $T_{i}(t) = T(x_{i}, t)$ is temperature
  • $w_{i}(t) = \rho_{w}(x_{i}, t)$ is the wood density
  • $c_{i}(t) = \rho_{c1}(x_{i}, t)$ is the char density
  • $a_{i}(t) = \rho_{c2}(x_{i}, t)$ is the ash density (if memory serves, this was another component in your kinetics from your previous question)
  • $i$ is an index for each grid point, so $i = 1, \ldots$ up to some number of grid points -- note that you'll need to handle boundary conditions appropriately
  • and you need to fill in the right-hand side with the appropriate values

If you can pose your equations in this way, then you can use built-in ODE solvers, which should enable you to attempt to solve these equations without having to explicitly form any matrices yourself (they will be assembled within any implicit ODE solver that calculates finite difference approximations of the Jacobian matrix). However, if you need to supply a Jacobian matrix, you can do so by differentiating the right-hand side of the ODEs above with respect to the dependent variables; this matrix is the one you're looking to form in your question.

  • $\begingroup$ I didn't think about using the ODE solvers in Matlab, which I guess I could try initially. What about other programming languages like C and Python, do they have ODE solvers? Eventually my model will need to run in C and Python, not just Matlab. That's why I'm trying to stay away from Matlab only functions. $\endgroup$
    – wigging
    Commented Sep 29, 2013 at 23:18
  • $\begingroup$ Yes, Python has several packages that solve ODEs. SciPy is the most famous (it does many other things too), but there are others. Similarly, many ODE solver packages exist in C. You could try SUNDIALS, PETSc, GSL, or write simple shims for packages written in Fortran (of which there are many). Some of the packages written in C have Python interfaces of varying quality, so it's possible to use the same solver across MATLAB, Python, and C, if you want. $\endgroup$ Commented Sep 30, 2013 at 0:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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