Implicit heat diffusion with kinetic reactions

Question

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?

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}

where:

• $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.