ZFK equation

Short description: Reaction–diffusion equation

ZFK equation, abbreviation for Zeldovich–Frank-Kamenetskii equation, is a reaction–diffusion equation that models premixed flame propagation. The equation is named after Yakov Zeldovich and David A. Frank-Kamenetskii who derived the equation in 1938 and is also known as the Nagumo equation. The equation is analogous to KPP equation except that is contains an exponential behaviour for the reaction term and it differs fundamentally from KPP equation with regards to the propagation velocity of the traveling wave. In non-dimensional form, the equation reads

[math]\displaystyle{ \frac{\partial \theta}{\partial t} = \frac{\partial^2\theta}{\partial x^2} + \omega(\theta) }[/math]

with a typical form for [math]\displaystyle{ \omega }[/math] given by

[math]\displaystyle{ \omega =\frac{\beta^2}{2} \theta(1-\theta) e^{-\beta(1-\theta)} }[/math]

where [math]\displaystyle{ \theta\in [0,1] }[/math] is the non-dimensional dependent variable (typically temperature) and [math]\displaystyle{ \beta }[/math] is the Zeldovich number. In the ZFK regime, [math]\displaystyle{ \beta\gg 1 }[/math]. The equation reduces to Fisher's equation for [math]\displaystyle{ \beta\ll 1 }[/math] and thus [math]\displaystyle{ \beta\ll 1 }[/math] corresponds to KPP regime. The minimum propagation velocity [math]\displaystyle{ U_{min} }[/math] (which is usually the long time asymptotic speed) of a traveling wave in the ZFK regime is given by

[math]\displaystyle{ U_{ZFK} \propto \sqrt{2\int_0^1\omega(\theta) d\theta} }[/math]

whereas in the KPP regime, it is given by

[math]\displaystyle{ U_{KPP} = 2 \sqrt{\left.\frac{d\omega}{d\theta}\right |_{\theta=0}}. }[/math]

Traveling wave solution

Numerical solution of ZFK equation

Similar to Fisher's equation, a traveling wave solution can be found for this problem. Suppose the wave to be traveling from right to left with a constant velocity [math]\displaystyle{ U }[/math], then in the coordinate attached to the wave, i.e., [math]\displaystyle{ z=x+Ut }[/math], the problem becomes steady. The ZFK equation reduces to

[math]\displaystyle{ U\frac{d\theta}{dz} = \frac{d^2\theta}{dz^2} + \frac{\beta^2}{2} \theta(1-\theta)e^{-\beta(1-\theta)} }[/math]

satisfying the boundary conditions [math]\displaystyle{ \theta(-\infty)=0 }[/math] and [math]\displaystyle{ \theta(\infty)=1 }[/math]. The boundary conditions are satisfied sufficiently smoothly so that the derivative [math]\displaystyle{ d\theta/dz }[/math] also vanishes as [math]\displaystyle{ z\rightarrow \pm\infty }[/math]. Since the equation is translationally invariant in the [math]\displaystyle{ z }[/math] direction, an additional condition, say for example [math]\displaystyle{ \theta(0)=1/2 }[/math], can be used to fix the location of the wave. The speed of the wave [math]\displaystyle{ U }[/math] is obtained as part of the solution, thus constituting a nonlinear eigenvalue problem. Numerical solution of the above equation, [math]\displaystyle{ \theta }[/math], the eigenvalue [math]\displaystyle{ U }[/math] and the corresponding reaction term [math]\displaystyle{ \omega }[/math] are shown in the figure, calculated for [math]\displaystyle{ \beta=15 }[/math].

Asymptotic solution

The ZFK regime as [math]\displaystyle{ \beta\rightarrow\infty }[/math] is formally analyzed using activation energy asymptotics. Since [math]\displaystyle{ \beta }[/math] is large, the term [math]\displaystyle{ e^{-\beta(1-\theta)} }[/math] will make the reaction term practically zero, however that term will be non-negligible if [math]\displaystyle{ 1-\theta \sim 1/\beta }[/math]. The reaction term will also vanish when [math]\displaystyle{ \theta=0 }[/math] and [math]\displaystyle{ \theta=1 }[/math]. Therefore, it is clear that [math]\displaystyle{ \omega }[/math] is negligible everywhere except in a thin layer close to the right boundary [math]\displaystyle{ \theta=1 }[/math]. Thus the problem is split into three regions, an inner diffusive-reactive region flanked on either side by two outer convective-diffusive regions.

Outer region

The problem for outer region is given by

[math]\displaystyle{ U\frac{d\theta}{dz} = \frac{d^2\theta}{dz^2}. }[/math]

The solution satisfying the condition [math]\displaystyle{ \theta(-\infty)=0 }[/math] is [math]\displaystyle{ \theta=e^{Uz} }[/math]. This solution is also made to satisfy [math]\displaystyle{ \theta(0)=1 }[/math] (an arbitrary choice) to fix the wave location somewhere in the domain because the problem is translationally invariant in the [math]\displaystyle{ z }[/math] direction. As [math]\displaystyle{ z\rightarrow 0^- }[/math], the outer solution behaves like [math]\displaystyle{ \theta=1+Uz + \cdots }[/math] which in turn implies [math]\displaystyle{ d\theta/dz=U + \cdots. }[/math]

The solution satisfying the condition [math]\displaystyle{ \theta(\infty)=1 }[/math] is [math]\displaystyle{ \theta=1 }[/math]. As [math]\displaystyle{ z\rightarrow 0^+ }[/math], the outer solution behaves like [math]\displaystyle{ \theta=1 }[/math] and thus [math]\displaystyle{ d\theta/dz=0 }[/math].

We can see that although [math]\displaystyle{ \theta }[/math] is continuous at [math]\displaystyle{ z=0 }[/math], [math]\displaystyle{ d\theta/dz }[/math] has a jump at [math]\displaystyle{ z=0 }[/math]. The transition between the derivatives is described by the inner region.

Inner region

In the inner region where [math]\displaystyle{ 1-\theta \sim 1/\beta }[/math], reaction term is no longer negligible. To investigate the inner layer structure, one introduces a stretched coordinate encompassing the point [math]\displaystyle{ z=0 }[/math] because that is where [math]\displaystyle{ \theta }[/math] is approaching unity according to the outer solution and a stretched dependent variable according to [math]\displaystyle{ \eta = \beta z, \, \Theta = \beta(1-\theta). }[/math] Substituting these variables into the governing equation and collecting only the leading order terms, we obtain

[math]\displaystyle{ 2\frac{d^2\Theta}{d\eta^2} = \Theta e^{-\Theta}. }[/math]

The boundary condition as [math]\displaystyle{ \eta\rightarrow -\infty }[/math] comes from the local behaviour of the outer solution obtained earlier, which when we write in terms of the inner zone coordinate becomes [math]\displaystyle{ \Theta \rightarrow -U\eta=+\infty }[/math] and [math]\displaystyle{ d\Theta/d\eta=-U }[/math]. Similarly, as [math]\displaystyle{ \eta\rightarrow\infty }[/math]. we find [math]\displaystyle{ \Theta=d\Theta/d\eta=0 }[/math]. The first integral of the above equation after imposing these boundary conditions becomes

[math]\displaystyle{ \begin{align} \left.\left(\frac{d\Theta}{d\eta}\right)^2\right |_{\Theta=\infty} - \left.\left(\frac{d\Theta}{d\eta}\right)^2\right |_{\Theta=0} &= \int_0^\infty \Theta e^{-\Theta}d\Theta\\ U^2 &= 1 \end{align} }[/math]

which implies [math]\displaystyle{ U=1 }[/math]. It is clear from the first integral, the wave speed square [math]\displaystyle{ U^2 }[/math] is proportional to integrated (with respect to [math]\displaystyle{ \theta }[/math]) value of [math]\displaystyle{ \omega }[/math] (of course, in the large [math]\displaystyle{ \beta }[/math] limit, only the inner zone contributes to this integral). The first integral after substituting [math]\displaystyle{ U=1 }[/math] is given by

[math]\displaystyle{ \frac{d\Theta}{d\eta}= - \sqrt{1-(\Theta+1)\exp(-\Theta)}. }[/math]

KPP–ZFK transition

In the KPP regime, [math]\displaystyle{ U_{min}=U_{KPP}. }[/math] For the reaction term used here, the KPP speed that is applicable for [math]\displaystyle{ \beta\ll 1 }[/math] is given by

[math]\displaystyle{ U_{KPP} = 2 \sqrt{\left.\frac{d\omega}{d\theta}\right |_{\theta=0}}= \sqrt 2 \beta e^{-\beta/2} }[/math]

whereas in the ZFK regime, as we have seen above [math]\displaystyle{ U_{ZFK}=1 }[/math]. Numerical integration of the equation for various values of [math]\displaystyle{ \beta }[/math] showed that there exists a critical value [math]\displaystyle{ \beta_*=1.64 }[/math] such that only for [math]\displaystyle{ \beta\leq \beta_* }[/math], [math]\displaystyle{ U_{min}=U_{KPP}. }[/math] For [math]\displaystyle{ \beta\geq \beta_* }[/math], [math]\displaystyle{ U_{min} }[/math] is greater than [math]\displaystyle{ U_{KPP} }[/math]. As [math]\displaystyle{ \beta\gg 1 }[/math], [math]\displaystyle{ U_{min} }[/math] approaches [math]\displaystyle{ U_{ZFK}=1 }[/math] thereby approaching the ZFK regime. The region between the KPP regime and the ZFK regime is called the KPP–ZFK transition zone.

The critical value depends on the reaction model, for example we obtain

[math]\displaystyle{ \beta_*=3.04 \quad \text{for}\quad \omega \propto (1-\theta) e^{-\beta(1-\theta)} }[/math]

[math]\displaystyle{ \beta_*=5.11 \quad \text{for}\quad \omega \propto (1-\theta)^2 e^{-\beta(1-\theta)}. }[/math]

Clavin–Liñán model

To predict the KPP–ZFK transition analytically, Paul Clavin and Amable Liñán proposed a simple piecewise linear model

[math]\displaystyle{ \omega(\theta)=\begin{cases} \theta \quad \text{if} \quad 0\leq \theta\leq 1-\epsilon,\\ h(1-\theta)/\epsilon^2 \quad \text{if} \quad 1-\epsilon\leq \theta\leq 1 \end{cases} }[/math]

where [math]\displaystyle{ h }[/math] and [math]\displaystyle{ \epsilon }[/math] are constants. The KPP velocity of the model is [math]\displaystyle{ U_{KPP}=2 }[/math], whereas the ZFK velocity is obtained as [math]\displaystyle{ U_{ZFK}=\sqrt h }[/math] in the double limit [math]\displaystyle{ \epsilon\rightarrow 0 }[/math] and [math]\displaystyle{ h\rightarrow\infty }[/math] that mimics a sharp increase in the reaction near [math]\displaystyle{ \theta=1 }[/math].

For this model there exists a critical value [math]\displaystyle{ h_*=1-\epsilon^2 }[/math] such that

[math]\displaystyle{ \begin{cases} h\lt h_*: &\quad U_{min}=U_{KPP},\\ h\gt h_*: &\quad U_{min}=\frac{h/(1-\epsilon)+1-\epsilon}{\sqrt{h/(1-\epsilon)-\epsilon}},\\ h\gg h_*: &\quad U_{min}\rightarrow U_{ZFK} \end{cases} }[/math]

See also

• Fisher's equation