Non-dimensionalization

Non-dimensional scales of time, temperature, length, and heat transfer may be defined as

${\displaystyle \tau ={\frac {t}{t_{e}}},\quad \theta ={\frac {\beta (T-T_{o})}{T_{o}}},\quad \eta ^{j}={\frac {r^{j}}{a}},\quad \delta ={\frac {t_{c}}{t_{e}}}}$

where

${\displaystyle t_{c}=\rho c_{v}a^{2}/\lambda }$	is the characteristic heat conduction time across the vessel;
${\displaystyle t_{f}=\left(Be^{-\beta }\right)^{\!-1}}$	is the characteristic fuel consumption time;
${\displaystyle t_{e}=\left(\beta \gamma Be^{-\beta }\right)^{\!-1}}$	is the characteristic explosion/ignition time;
${\displaystyle a}$	is the characteristic distance, e.g., vessel radius;
${\displaystyle \beta =E/(R\,T_{o})}$	is the non-dimensional activation energy;
${\displaystyle \gamma =(qY_{Fo})/(c_{v}T_{o})}$	is the heat-release parameter;
${\displaystyle \delta }$	is the Damköhler number;
${\displaystyle r}$	is the spatial coordinate with origin at the center;
${\displaystyle j=0}$	for planar slab;
${\displaystyle j=1}$	for cylindrical vessel;
${\displaystyle j=2}$	for spherical vessel.

Note

In a typical combustion process, ${\displaystyle \gamma \approx 6{-}8,\ \beta \approx 30{-}100}$ so that ${\displaystyle \beta \gamma \gg 1}$.

Therefore, ${\displaystyle t_{f}=\beta \gamma t_{e}\gg 1}$. That is, fuel consumption time is much longer than ignition time, so fuel consumption is essentially negligible in the study of ignition.

This is why the fuel concentration is assumed to remain the initial fuel concentration ${\displaystyle Y_{Fo}}$.

Substituting the non-dimensional variables in the energy equation from the introduction

${\displaystyle {\frac {\partial \theta }{\partial \tau }}={\frac {1}{\delta }}{\frac {1}{\eta ^{j}}}{\frac {\partial }{\partial \eta }}\left(\eta ^{j}{\frac {\partial \theta }{\partial \eta }}\right)+e^{\theta /(1+\theta /\beta )}}$

Since ${\displaystyle \beta \gg 1}$, the exponential term can be linearized ${\displaystyle e^{\theta /(1+\theta /\beta )}\approx e^{\theta }}$, hence

${\displaystyle {\frac {\partial \theta }{\partial \tau }}={\frac {1}{\delta }}{\frac {1}{\eta ^{j}}}{\frac {\partial }{\partial \eta }}\left(\eta ^{j}{\frac {\partial \theta }{\partial \eta }}\right)+e^{\theta }}$

At ${\displaystyle \tau =0}$, we have ${\displaystyle \theta (\eta ,0)=0}$ and for ${\displaystyle \tau >0}$, ${\displaystyle \theta }$ needs to satisfy ${\displaystyle \theta (1,\tau )=0}$ and ${\displaystyle \partial \theta /\partial \eta |_{\eta =0}=0.}$

Steady-state regime

When ${\displaystyle 0<\delta <e^{-1}}$, the linear term eventually dominates and the system is able to reach a steady state as ${\displaystyle \tau \rightarrow \infty }$. At steady state (${\displaystyle d\theta /d\tau =0}$), the balance is given by the equation

${\displaystyle \delta e^{\theta }=\theta ,\quad \Rightarrow \quad \theta =-W(-\delta )}$

where ${\displaystyle W}$ represents the Lambert W function. From the properties of Lambert W function, it is easy to see that the steady state temperature provided by the above equation exists only when ${\displaystyle \delta \leq \delta _{c}=1/e}$, where ${\displaystyle \delta _{c}}$ is called as Frank-Kamenetskii parameter as a critical point where the system bifurcates from the existence of steady state to explosive state at large times.

Explosive regime

For ${\displaystyle \delta _{c}<\delta <\infty }$, the system explodes since the exponential term dominates as time proceeds. We do not need to wait for a long time for ${\displaystyle \theta }$ to blow up. Because of the exponential forcing, ${\displaystyle \theta \rightarrow \infty }$ at a finite value of ${\displaystyle \tau }$. This time is interpreted as the ignition time or induction time of the system. When ${\displaystyle \delta \gg \delta _{c}}$, the heat conduction term ${\displaystyle -\theta /\delta }$ can be neglected in which case the problem admits an explicit solution,

${\displaystyle {\frac {d\theta }{d\tau }}=e^{\theta },\quad \Rightarrow \quad \theta =\ln \left({\frac {1}{1-\tau }}\right)}$

At time ${\displaystyle \tau =1}$, the system explodes. This time is also referred to as the adiabatic induction period since the heat conduction term ${\displaystyle -\theta /\delta }$ is neglected.

In the near-critical condition, i.e., when ${\displaystyle \delta -\delta _{c}\rightarrow 0^{+}}$, the system takes very long time to explode. The analysis for this limit was first carried out by Frank-Kamenetskii.,^[10] although proper asymptotics were carried out only later by D. R. Kassoy and Amable Liñán^[11] including reactant consumption because reactant consumption is not negligible when ${\displaystyle \tau \sim \beta \gamma }$. A simplified analysis without reactant consumption is presented here. Let us define a small parameter ${\displaystyle \epsilon \rightarrow 0^{+}}$ such that ${\displaystyle \delta =\delta _{c}(1+\epsilon )}$. For this case, the time evolution of ${\displaystyle \theta }$ is as follows: first it increases to steady-state temperature value corresponding to ${\displaystyle \delta =\delta _{c}}$, which is given by ${\displaystyle \theta =-W(-\delta _{c})=1}$ at times of order ${\displaystyle \tau \sim O(1)}$, then it stays very close to this steady-state value for a long time before eventually exploding at a long time. The quantity of interest is the long-time estimate for the explosion. To find out the estimate, introduce the transformations ${\displaystyle \zeta ={\sqrt {\epsilon }}\tau }$ and ${\displaystyle \theta =1+{\sqrt {\epsilon }}\psi (\zeta )-\epsilon }$ that is appropriate for the region where ${\displaystyle \theta }$ stays close to ${\displaystyle 1}$ into the governing equation and collect only the leading-order terms to find out

${\displaystyle \delta _{c}{\frac {d\psi }{d\zeta }}=1+{\frac {\psi ^{2}}{2}},\quad \psi (0)\rightarrow -(1-\theta )/{\sqrt {\epsilon }}=-\infty }$

where the boundary condition is derived by matching with the initial region wherein ${\displaystyle \tau ,\theta \sim 1}$. The solution to the above-mentioned problem is given by

${\displaystyle \psi =-{\sqrt {2}}\cot \left({\frac {\zeta }{{\sqrt {2}}\delta _{c}}}\right)}$

which immediately reveals that ${\displaystyle \psi \rightarrow \infty }$ when ${\displaystyle \zeta ={\sqrt {2}}\pi \delta _{c}.}$ Writing this condition in terms of ${\displaystyle \tau }$, the explosion time in the near-critical condition is found to be

${\displaystyle \tau ={\sqrt {\frac {2\pi ^{2}\delta _{c}^{3}}{\delta -\delta _{c}}}}}$

which implies that the ignition time ${\displaystyle \tau \rightarrow \infty }$ as ${\displaystyle \delta -\delta _{c}\rightarrow 0^{+}}$ with a square-root singularity.

Planar vessel

Frank-Kamenetskii explosion for planar vessel

For planar vessel, there is an exact solution. Here ${\displaystyle j=0}$, then

${\displaystyle {\frac {d^{2}\theta }{d\eta ^{2}}}=-\delta e^{\theta }}$

If the transformations ${\displaystyle \Theta =\theta _{m}-\theta }$ and ${\displaystyle \xi ^{2}=\delta e^{\theta _{m}}\eta ^{2}}$, where ${\displaystyle \theta _{m}}$ is the maximum temperature which occurs at ${\displaystyle \eta =0}$ due to symmetry, are introduced

${\displaystyle {\frac {d^{2}\Theta }{d\xi ^{2}}}=e^{-\Theta },\quad \Theta (0)=0,\quad {\frac {d\Theta }{d\xi }}_{\xi =0}=0}$

Integrating once and using the second boundary condition, the equation becomes

${\displaystyle {\frac {d\Theta }{d\xi }}={\sqrt {2(1-e^{-\Theta })}}}$

and integrating again

${\displaystyle e^{(\theta _{m}-\theta )/2}=\cosh \left(\eta {\sqrt {\frac {\delta e^{\theta _{m}}}{2}}}\right)}$

The above equation is the exact solution, but ${\displaystyle \theta _{m}}$ maximum temperature is unknown, but we have not used the boundary condition of the wall yet. Thus using the wall boundary condition ${\displaystyle \theta =0}$ at ${\displaystyle \eta =1}$, the maximum temperature is obtained from an implicit expression,

${\displaystyle e^{\theta _{m}/2}=\cosh {\sqrt {\frac {\delta e^{\theta _{m}}}{2}}}\quad {\text{or}}\quad \delta =2e^{-\theta _{m}}\left(\operatorname {arcosh} e^{\theta _{m}/2}\right)^{\!2}}$

Critical ${\displaystyle \delta _{c}}$ is obtained by finding the maximum point of the equation (see figure), i.e., ${\displaystyle d\delta /d\theta _{m}=0}$ at ${\displaystyle \delta _{c}}$.

${\displaystyle {\frac {d\delta }{d\theta _{m}}}=0,\quad \Rightarrow \quad e^{\theta _{m,c}/2}-{\sqrt {e^{\theta _{m,c}}-1}}\operatorname {arcosh} e^{\theta _{m,c}/2}=0}$

${\displaystyle \theta _{m,c}=1.1868,\quad \Rightarrow \quad \delta _{c}=0.8785}$

So the critical Frank-Kamentskii parameter is ${\displaystyle \delta _{c}=0.8785}$. The system has no steady state (or explodes) for ${\displaystyle \delta >\delta _{c}=0.8785}$ and for ${\displaystyle \delta <\delta _{c}=0.8785}$, the system goes to a steady state with very slow reaction.

Cylindrical vessel

Frank-Kamenetskii explosion for cylindrical vessel

For cylindrical vessel, there is an exact solution. Though Frank-Kamentskii used numerical integration assuming there is no explicit solution, Paul L. Chambré provided an exact solution in 1952.^[14] H. Lemke also solved provided a solution in a somewhat different form in 1913.^[15] Here ${\displaystyle j=1}$, then

${\displaystyle {\frac {1}{\eta }}{\frac {d}{d\eta }}\left(\eta {\frac {d\theta }{d\eta }}\right)=-\delta e^{\theta }}$

If the transformations ${\displaystyle \omega =\eta d\theta /d\eta }$ and ${\displaystyle \chi =e^{\theta }\eta ^{2}}$ are introduced

${\displaystyle {\frac {d\omega }{d\chi }}=-{\frac {\delta }{2+\omega }},\quad \omega (0)=0}$

The general solution is ${\displaystyle \omega ^{2}+4\omega +C=-2\delta \chi }$. But ${\displaystyle C=0}$ from the symmetry condition at the centre. Writing back in original variable, the equation reads,

${\displaystyle \eta ^{2}\left({\frac {d\theta }{d\eta }}\right)^{2}+4\eta {\frac {d\theta }{d\eta }}=-2\delta \eta ^{2}e^{\theta }}$

But the original equation multiplied by ${\displaystyle 2\eta ^{2}}$ is

${\displaystyle 2\eta ^{2}{\frac {d^{2}\theta }{d\eta ^{2}}}+2\eta {\frac {d\theta }{d\eta }}=-2\delta \eta ^{2}e^{\theta }}$

Now subtracting the last two equation from one another leads to

${\displaystyle {\frac {d^{2}\theta }{d\eta ^{2}}}-{\frac {1}{\eta }}{\frac {d\theta }{d\eta }}-{\frac {1}{2}}\left({\frac {d\theta }{d\eta }}\right)^{2}=0}$

This equation is easy to solve because it involves only the derivatives, so letting ${\displaystyle g(\eta )=d\theta /d\eta }$ transforms the equation

${\displaystyle {\frac {dg}{d\eta }}-{\frac {g}{\eta }}-{\frac {g^{2}}{2}}=0}$

This is a Bernoulli differential equation of order ${\displaystyle 2}$, a type of Riccati equation. The solution is

${\displaystyle g(\eta )={\frac {d\theta }{d\eta }}=-{\frac {4\eta }{B'+\eta ^{2}}}}$

Integrating once again, we have ${\displaystyle \theta =A-2\ln(B\eta ^{2}+1)}$ where ${\displaystyle B=1/B'}$. We have used already one boundary condition, there is one more boundary condition left, but with two constants ${\displaystyle A,\ B}$. It turns out ${\displaystyle A}$ and ${\displaystyle B}$ are related to each other, which is obtained by substituting the above solution into the starting equation we arrive at ${\displaystyle A=\ln(8B/\delta )}$. Therefore, the solution is

${\displaystyle \theta =\ln {\frac {8B/\delta }{(B\eta ^{2}+1)^{2}}}}$

Now if we use the other boundary condition ${\displaystyle \theta (1)=0}$, we get an equation for ${\displaystyle B}$ as ${\displaystyle \delta (B+1)^{2}-8B=0}$. The maximum value of ${\displaystyle \delta }$ for which solution is possible is when ${\displaystyle B=1}$, so the critical Frank-Kamentskii parameter is ${\displaystyle \delta _{c}=2}$. The system has no steady state( or explodes) for ${\displaystyle \delta >\delta _{c}=2}$ and for ${\displaystyle \delta <\delta _{c}=2}$, the system goes to a steady state with very slow reaction. The maximum temperature ${\displaystyle \theta _{m}}$ occurs at ${\displaystyle \eta =0}$

${\displaystyle \theta _{m}=\ln {\frac {8B}{\delta }}\quad {\text{or}}\quad \delta =8Be^{-\theta _{m}}}$

For each value of ${\displaystyle \delta }$, we have two values of ${\displaystyle \theta _{m}}$ since ${\displaystyle B}$ is multi-valued. The maximum critical temperature is ${\displaystyle \theta _{m,c}=\ln 4}$.

Spherical vessel

Frank-Kamenetskii spherical vessel

For spherical vessel, there is no known explicit solution, so Frank-Kamenetskii used numerical methods to find the critical value. Here ${\displaystyle j=2}$, then

${\displaystyle {\frac {1}{\eta ^{2}}}{\frac {d}{d\eta }}\left(\eta ^{2}{\frac {d\theta }{d\eta }}\right)=-\delta e^{\theta }}$

If the transformations ${\displaystyle \Theta =\theta _{m}-\theta }$ and ${\displaystyle \xi ^{2}=\delta e^{\theta _{m}}\eta ^{2}}$, where ${\displaystyle \theta _{m}}$ is the maximum temperature which occurs at ${\displaystyle \eta =0}$ due to symmetry, are introduced

${\displaystyle {\frac {1}{\xi ^{2}}}{\frac {d}{d\xi }}\left(\xi ^{2}{\frac {d\Theta }{d\xi }}\right)=e^{-\Theta },\quad \Theta (0)=0,\quad {\frac {d\Theta }{d\xi }}_{\xi =0}=0}$

The above equation is nothing but Emden–Chandrasekhar equation,^[16] which appears in astrophysics describing isothermal gas sphere. Unlike planar and cylindrical case, the spherical vessel has infinitely many solutions for ${\displaystyle \delta <\delta _{c}}$ oscillating about the point ${\displaystyle \delta =2}$,^[17] instead of just two solutions, which was shown by Israel Gelfand.^[18] The lowest branch will be chosen to explain explosive behavior.

From numerical solution, it is found that the critical Frank-Kamenetskii parameter is ${\displaystyle \delta _{c}=3.3220}$. The system has no steady state( or explodes) for ${\displaystyle \delta >\delta _{c}=3.3220}$ and for ${\displaystyle \delta <\delta _{c}=3.3220}$, the system goes to a steady state with very slow reaction. The maximum temperature ${\displaystyle \theta _{m}}$ occurs at ${\displaystyle \eta =0}$ and maximum critical temperature is ${\displaystyle \theta _{m,c}=1.6079}$.