数理解析研究所講究録
1191 巻 2001 年 43-50

43Numerical simulation of helical waves arising in
self-propagating high-temperature syntheses
Tsutomu Ikeda (Ryukoku Univ., tsutomu@rins.ryukoku.ac.jp)
Masaharu Nagayama (Kyoto Univ., nagayama@kurims.kyoto-u.ac.jp)
池田勉 (龍谷大学)

長山雅晴 (京都大学)

Self-propagating High-temperature Syntheses (SHS)
Self-propagating High-temperature Syntheses (SHS) recently attract a great deal ofattention as a synthetic method creating refractory ceramics, intermetalic compounds,
composites and so on ([9] for instance). Among over 100 kinds of SHS, a typical example
is the synthesis of titanium carbide
from titanium Ti and carbon C. In the process,we ignite the top end of cylindrical sample formed by compressing powdery raw materials,
say Ti and , and the temperature there goes up to over 3,000 by a great quantity ofheat emitted along with chemical reactions between raw materials in the sample. Thus a
combustion disk appears and the synthesis spreads to the neighborhood warmed by the heatdiffusion effect. Once the top end of compact is set fire to, then a reaction front propagatesalong the compact with almost no energy $\sup$.ply and the synthesis in the compact of 4 cm
long completes in a few seconds.
When a combustion wave keeps its profile and propagates at a constant velocity, one
can create a very-high-quality uniform product by the SHS. In other words, combustion
waves of steady-state mode bring high-quality products. When experimental conditions
are changed, however, the combustion wave of steady-state mode may lose its stability andgives place some non-uniform combustion wave, which deteriorates quality. Actually, it
is reported that a pulsating combustion wave emerges through the Hopf bifurcation of a
combustion wave of steady-state mode ([9]). We also observe a wave that propagates in
the form of spiral encircling the cylindrical sample with several reaction spots, which iscalled a helical combustion wave in the present paper. We note that similar helical waves
are reported in autocatalytic reactions and propagation fronts of polymerizations (Pojman
et al. [10] .
In the present paper, by using mathematical models for the SHS process, we try to
clarify the existing condition of stable helical wave and transition processes of wave patterns. We have studied numericallyfrom steady-state mode to helical mode via pulsating
the behaviors of propagating waves in cylindrical domains and have obtained the following
results:(a) Set physical parameters so that a pulsating wave exists stably for the one-dimensional
problem. Then, the planar pulsating wave is still stable in the cylindrical domain when
the radius is small while a helical wave takes the place of the pulsating wave when the theradius becomes larger.
(b) When physical parameters are fixed so that a combustion wave of steady-state mode isstable in the one-dimensional problem, the planar wave of steady-state mode is stable also$\mathrm{T}\mathrm{i}\mathrm{C}$$\mathrm{C}$$\mathrm{K}$$)$$\mathrm{o}\mathrm{n}\mathrm{e}\dot{\mathrm{s}}$44in the cylindrical domain.

A dimensionless mathematical model for SHS
We begin with a mathematical model for a two-component spatially uniform synthesis$m \mathrm{R}_{1}+n\mathrm{R}_{2}\frac{k(T)}{},$$\mathrm{P}$

(1)

,

where
and
denote reactants, a product and $m$ and reaction orders. The chemical
reaction rate $k(T)$ is a function of temperature governed by the Arrhenius law, and it isexpressed as$\mathrm{P}$$\mathrm{R}_{2}$$\mathrm{R}_{1}$

$n$

$T$

$k(T)=k\mathrm{o}\emptyset(\tau)$,$\phi(T)=\{$$\exp(-E/RT)$$T\geq T_{i_{\mathit{9}}}$$0$$\tau<\tau_{ig}$,(2)by the use of a pre-exponential constant , an effective activation energy , the universal
gas constant $R=$ 8.3135
and an ignition temperature. This form of $k(T)$means that no reaction occurs when is less thanand the synthesis starts when goes$E$

$k_{0}$$\mathrm{J}/\mathrm{m}\mathrm{o}1\cdot \mathrm{K}$$T_{ig}$

$T$

up to
. Under the assumption
concentration of reactants
andby$T_{ig}$$\mathrm{R}_{1}$$C_{p}\rho\dot{T}$$T$$T_{ig}$, whereanddenote the initial, respectively, the model system of equations is given$\frac{1}{m}[\mathrm{R}_{1}]=\frac{1}{m}[\mathrm{R}_{2}]$$\mathrm{R}_{2}$$[\mathrm{R}_{1}]$$[\mathrm{R}_{2}]$,$=$$Qk_{0^{m^{m}}}n^{n_{C\emptyset}}m+n(\tau)$$\dot{c}=$$-k_{0^{m^{mn_{C^{m+n}}}}}n\emptyset(\tau)$(3)$t>0$ ,$\{$.

where the dot means the differentiation by time , and the unknown variables are the$t$temperature $T=T(t)$ and(the concentration of ) at time . Here,denotesthe specific heat capacity, the density of compact and the reaction heat.
In general, bothand depend on the temperature, however, we assume them to be
paper by taking their suitable average. Thus, integrating the bothconstant in the$C_{p}\rho
T=-Q\dot{c}$, we can define the burned temperature
hand side ofby$c(t)= \frac{1}{m}$$t$$\mathrm{R}_{1}$$C_{p}$$Q$$\rho$$Q$$C_{p}$$\mathrm{p}\mathrm{r}.\mathrm{e}\mathrm{s}\mathrm{e}\mathrm{n}\mathrm{t}$$T_{b}$(4)$T_{b}= \tau_{0+}\frac{Q}{C_{p}\rho}C0$

where $T_{0}=T(0)$ and $c_{0}=c(0)$ are the initial temperature and concentration, respectively.
We note that the solution $(T(t), c(t))$ of (3) tends toasif.
Considering the spatial distribution $T(t, x)$ and $c(t, x)$ of temperature and reactant concentration and introducing the heat diffusion effect, we arrive at the following mathematicalmodel for SHS:
$tarrow\infty$$(T_{b}, 0)$$C_{p} \rho\frac{\partial T}{\partial t}$$=$$C_{p}\rho\alpha\triangle T+Qk_{0}m^{m}nc\emptyset nm+n(T)$

,

in

$\{$

$\frac{\partial c}{\partial t}$$=$$-k_{0^{mn}}mnC^{m+}\phi n(T)$,$T_{0}\geq T_{ig}$$\Omega$, $t>0$ ,

(5)

45

where
is a diffusion constant anddenotes the cylindrical domain occupied by the
sample compact. Introducing dimensionless parameters and variables defined by$\alpha$$\Omega$$K0=k_{00}m^{m}nC^{m}n+n-1$ ,$\tau=K_{0}t$,,$y=\sqrt{\frac{K_{0}}{\alpha}}x$

(6)

$\{$

$u( \tau, y)=\frac{T(\frac{\tau}{K_{0}},\sqrt{\frac{\alpha}{K_{0}}}y)-\tau 0}{T_{b}-T_{0}}$,$\eta(\tau, y)=1-\frac{1}{c_{0}}C(\frac{\tau}{K_{0}}, \sqrt{\frac{\alpha}{K_{0}}}y)$and rewriting
and $\sqrt{K_{0}/\alpha}\Omega=\{\sqrt{K_{0}/\alpha}x;x\in\Omega\}$ respectively as
derive the following non-dimensional system from (5):$\tau,$$y$$=$$\frac{\partial u}{\partial t}$

in$=$$(.1-\eta)N\Phi(u)$

$x$

and

$\Omega$, we,$\triangle u+(1-\eta)^{N}\Phi(u)$$\{$$\frac{\partial\eta}{\partial t}$$t,$$\Omega$, $t>0$ .

(7),

The new variables and are called non-dimensional temperature and a synthesis rate,
respectively. In (7), $N=m+n$ andis expressed as$u$$\eta$$\Phi(u)$$\exp(-\frac{e}{u+u_{0}})$$u\geq u_{i_{\mathit{9}}}$$\Phi(u)=\{$$0$$u<u_{ig}$

,

(8)

,

by using non-dimensional parameters$e= \frac{E}{R(T_{b}-T_{0})}$,$u_{0}= \frac{T_{0}}{T_{b}-T_{0}}$When the radius of cylindrical domain
dimensional problem:$\Omega$,tends to$\frac{\partial u}{\partial t}=\frac{\partial^{2}u}{\partial z^{2}}+(1-\eta)^{N}\Phi(u)$,$u_{i_{\mathit{9}}}= \frac{T_{i_{\mathit{9}^{-}}}\tau_{0}}{T_{b}-T_{0}}$

$0$

.

(9)

, (7) may be reduced to the one-$\frac{\partial\eta}{\partial t}=(1-\eta)N\Phi(u)$

.

(10)

A solution $\{u(t, z), \eta(t, z)\}$ of (10) is called a travelling wave when it is of the form$u(t, z)=\tilde{u}(z-ct)$ and $\eta(t, z)=\tilde{\eta}(z-ct)$ withconstant . A travelling wave corresponds to a combustion wave of steady-state mode, that is, it keeps the profileand propagates at the constant velocity .
$c$$\{\tilde{u}(\cdot),\tilde{\eta}(\cdot)\}$$c$

Review of mathematical results
We here briefly review mathematical results related to the gasless combustion problem(7). The existence of travelling wave of (10) has been shown by Logak and Loubeau[6] without stability analysis, and Bonnet and Logak [2] has implied the instability oftravelling wave for high effective activation energy in the one-dimensional problem (10).
$e$

46

There are many numerical studies on the one-dimensional problem, where is adopted asa paramete..r. The following transitions of propagating patterns have been reported as isincreased: appearance of pulsating waves through the Hopf bifurcation of travelling waves,chaotic behavior after period doubling bifurcations of pulsating waves ([1]), intermittent
propagation phenomena ([4]) and so on.
Ivleva et al. [5] have proposed a two-dimensional system of equations on the cylindricalsurface, and they have succeeded in numerical computation of two-dimensional helical waves
based on their system. On the other hand, Matkowsky and Sivashinsky [7] have shown
the destabilization of travelling wave via Hopf bifurcation by using their reduced equation
obtained from (5) in the limit where the effective activation energy $E$ tends to the infinity.
By the use of the reduced equation, Sivashinsky [11] has implied also the appearance ofhelical wave through a Hopf-type bifurcation of planar propagating wave.
As stated in the above, the existence of helical wave has been shown only for the twodimensional problem and the reduced equation. Different from [5] nor [11], we deal with thethree-dimensional full system of equations (7) in the present paper. Through our approach,
we will obtain the following new insights:(1) propagation patterns in the interior of cylindrical domain,(2) difference in propagation patterns between the two- and three-dimensional problems,(3) dependence of propagation patterns on the effective activation energy $E$ .
$e$

$e$

Numerical simulation
The reaction order $N$ should equal $m+n$ doubtless if the chemical reaction formula(1) exactly describes the SHS process, however, (1) is a highly idealized and simplified
for instance, there can occur the dissolution of carbon
one. In the actual synthesis oflayer. In order tointo a titanium melt and the diffusion of carbon through a solid$N$
introduce effects ignored in (1) into our mathematical model, we leave as a free parameter
independent of $m$ and .
In our numerical simulation of (7) (8) we decide values of parameters according to
combustion synthesis ([3] and [12]):perimental data of the$\mathrm{T}\mathrm{i}\mathrm{C}$$\mathrm{T}\mathrm{i}\mathrm{C}$

$n$

-

ex.

$\mathrm{T}\mathrm{i}\mathrm{C}$$m=1$ ,

$n=1$ ,

$N=3.0$ ,$T_{0}=298\mathrm{K}$,$T_{i_{\mathit{9}}}=1924\mathrm{K}$

,

(11)

$\{$

$C_{p}=53.78\mathrm{J}/\mathrm{m}\mathrm{o}1\cdot \mathrm{K}$,$Q_{M}=Qc_{0}/\rho=184.5\cross 10^{3}\mathrm{J}/\mathrm{m}\mathrm{o}1$.

are the averaged one from 298 to 3200 K. The value
andThe above values of$E$, and it depends on the particle sizes
energyisof orderof effective activation$E$
as a controllable parameter in our numericalof reactants. For this reason we leave
simulation.
$C_{p}$$Q_{M}$$\mathrm{K}$$10^{5}\mathrm{J}/\mathrm{m}\mathrm{o}1$One-dimensional patterns
In numerical simulations of the one-dimensional problem (10), as reported in previous
works, we also have observed stable travelling waves for small values of effective activationenergy , stable pulsating waves arising via the Hopf bifurcation of travelling waves for
$E$

47

medium $E$ , pulsating waves after period doubling bifurcation for large $E$ and chaotic combustion waves for much larger $E$ . We note that as $E$ increases the reaction stagnates andthe average velocity of propagation becomes slower.

Three-dimensional patterns
We consider (7) in the cylindrical domain $\Omega=\{(r, \theta, \mathcal{Z});0\leq r<R_{0},0\leq\theta<2\pi,$ $0<$$z<L_{z}\}$ subject to the boundary condition:$\frac{\partial u}{\partial z}=0$

for $z=0,$ $z=L_{z}$ and $t>0$ ,$\frac{\partial u}{\partial r}=0$

for.

$r=R_{0}$and $t>0$ .

(12)When the value of effective activation energy $E$ is set so that the one-dimensional problem(10) has a stable travelling wave, the planar travelling wave is also stable in the cylindrical
domain with various radiuses
in all of our numerical simulations.
$E$
We turn to the case where is selected so that (10) has a stable pulsating wave. The
planar pulsating wave is still stable in the cylindrical domain if
is small, however, it loses
its stability whenis large enough and helical combustions become dominant for large. Every helical wave has one or more reaction spots where the temperature is very high.
In some cases reaction spots seem to move independently and there occur annihilationof reaction spots by collision and creation of reaction spots (irregular helical waves). In
other cases several reaction spots seem to rotate on the cylindrical surface in the same
direction (regular helical wave). For the same physical and chemical parameters, there
coexist irregular patterns and regular patterns with various numbers of reaction spots.
These differences in patterns are caused only by differences in the initial distribution oftemperature. The rotation direction also depends on the initial distribution. Among
various helical patterns, an irregular one may be most popular in the sense that one has
to carefully set the initial distribution in order to get a regular helicalwhile oneneeds no special care for obtaining an irregular one.
In the isothermal surface of irregular wave, there appear elevation and subsidence erratically as time marches on. An isothermal surface of regular helical wave has some armsof which number is the same as that of reaction spots. The isothermal surface with arms
helically rotates down in the cylindrical domain as time passes on. Observing the distribution of temperature on horizontal cross section moving with the corresponding helicalwave, we know that a pulsating motion is dominant near the center while a helical rotation
is dominant near the cylindrical surface.
$R_{0}$$R_{0}$$R_{0}$$R_{0}$$\mathrm{P}^{\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}}\mathrm{r}\mathrm{n}\backslash$

Vestiges of helical combustion
The system (7) subject to the boundary condition (12) has helical waves as stated in
the above. However, because (12) means that there is no heat loss through the synthesis
process, $u(t, x)$ andtend to 1 as $tarrow\infty$ for any $x\in\Omega$ under suitable initialconditions and there is no vestige of helical combustion left at the completion of synthesis.
In the actual synthesis, it is hard to keep out heat loss from the cylindrical surface causedby convection and radiation. With a view to checking an effect of heat loss we replace the$\eta(t, x)$

48

second condition in (12) with the following simple heat loss condition by convection$\frac{\partial u}{\partial r}=-\lambda u$(13)

for $r=R_{0}$ and $t>0$ ,where is a positive constant. Just like the case of no heat loss boundary condition, thetemperature distribution becomes spatially uniform as $tarrow\infty$ . Nevertheless, vestiges ofhelical combustion are remained in the distribution of synthesis rate . We thus succeed in
the reproduction of vestiges of helical combustion similar to those observed in laboratory.
$\lambda$$\eta$

Differences between two- and three-dimensional patterns
We have carried out numerical simulations of the following plain two-dimensional problem$\frac{\partial u}{\partial t}=\frac{\partial^{2}u}{\partial\theta^{2}}+\frac{\partial^{2}u}{\partial z^{2}}+(1-\eta)^{N}\Phi(u)$for,$\frac{\partial\eta}{\partial t}--(1-\eta)^{N}\Phi(u)$$z\in(0, L_{z})$$\theta\in(0, L_{\theta}),$(14)

and $t>0$

subject to the periodic boundary condition in the -direction and the no-flux boundary
condition at.. $z=0$ and $z=L_{z}$ . The system (14) has helical waves in the same way as thethree-dimensional problem, that is, helical waves become the most popular type of solutions
is large enough. Differentof (14) when $E$ is set so that (10) has a pulsating wave andfrom the three-dimensional case, the two-dimensional system prefers regular helical waves
in the sense that solutions of (14) tend to a regular helical wave as time goes on for a verylarge family of initial conditions.
$\theta$$L_{\theta}$

Discussion.

One of our purpose in the present paper is to make clear the organizing center ofhelically propagating waves. We propose a conjecture that the organizing center is theexistence of stable pulsating wave of the corresponding one-dimensional problem, althougha mathematically rigorous proof has not been obtained. This conjecture is true for themodel (7) as described in the above. With the view of confirming the conjecture, we
consider an autocatalytic reaction model. The one-dimensional model is written as$\frac{\partial u}{\partial t}=\frac{\partial^{2}u}{\partial z^{2}}+vf(u)$,$\frac{\partial v}{\partial t}=d\frac{\partial^{2}v}{\partial z^{2}}-vf(u)$

,

(15)

where $f(u)=u^{m}(m\geq 1)$ and $d\geq 0$ is a constant, and Metcalf et al. [8] have proved that(15) has travelling waves in the case of $d=0$ . Metcalf et al. have also studied numericallythe behavior of solutions of (15) for various $m$ and , and it has been reported that (15)
has a stable travelling wave for some $m$ and while a travelling wave loses its stability
and a pulsating wave emerges for some other $m$ and . We have carried out numericalsimulation of the two-dimensional autocatalytic reaction model:
$d$

$d$

$d$

$\frac{\partial u}{\partial t}=\frac{\partial^{2}u}{\partial\theta^{2}}+\frac{\partial^{2}u}{\partial z^{2}}+vf(u)$for,$\frac{\partial v}{\partial t}=d(\frac{\partial^{2}v}{\partial\theta^{2}}+\frac{\partial^{2}v}{\partial z^{2}})-vf(u)$$\theta\in(0, L_{\theta}),$(16)$z\in(0, L_{z})$and $t>0$49

subject to the periodic boundary condition in the -direction and the no-flux boundary
condition at $z=0$ and $z=L_{z}$ , and have found the following properties supporting theabove mentioned conjecture:(a) When (15) has a stable pulsating wave, the planar pulsating wave is still stable in
the two-dimensional domain
ifis small while a helical wave becomesdominant for large .
(b) When (15) has a stable travelling wave, the planar travellingwave.is also stable in thetwo-dimensional domain with all the values of
we tried.
$\theta$$(0, L_{\theta})\cross(0, L_{z})$$L_{\theta}$$L_{\theta}$$L_{\theta}$

Acknowledgments
The present work was supported by the High-Tech Research Center (HRC) of Ryukoku
University. The present research fruits were obtained when the second author studied in the
HRC as a post-doctoral fellow. The authors thank the-Aid for Scientific Research(Grant No. 12440032, 12740062 and 12304006).
$\mathrm{G}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{t}_{-}\mathrm{i}\mathrm{n}$References[1] A. Bayliss and B. J. Matkowsky “Two routes to chaos in condensed phase combustion”,
SIAM J. Appl. Math. 50 (2), 437-459 (1990).

[2] A. Bonnet and E. Logak “Instability of travelling waves in solid combustion for highactivation energy”, to appear.
[3] S. D. Dunmead, D. W. Readey and C. E. Semler, “Kinetics of Combustion Synthesisin the Ti-C and
Systems”, J. Am. Ceram. Soc. 72 (12), 2318-2324 (1989).
$\mathrm{T}\mathrm{i}_{-}\mathrm{c}_{-}\mathrm{N}\mathrm{i}$[4] T. Ishiwata, M. Nagayama, T. Ikeda and M. Ohyanagi, “Numerical research of propagating form of combustion waves arising in a model for Self-propagating High temperature synthesis”, in preparation.
[5] T. P. Ivleva, A. G. Merzhanov and K. G. Shkadinskii, “Mathematical model of spin
combustion”, Sov. Phys. Dokl. 23 (4), 255-254 (1978).
[6] E. Logak and V. Loubeau, “Travelling waves solutions to a condensed phase combustion model”, Asymptotic Analysis 12, 259-294 (1996).

[7] B. J. Matkowsky and G. I. Sivashinsky, “Propagation of a pulsating reaction front insolid fuel combustion”, SIAM J. Appl. Math., 35, 465-478 (1978).
[8] M. J. Metcalf, J. H. Merkin and S. K. Scott, “Oscillating wave fronts in isothermal
chemical systems with arbitrary powers of autocatalysis”, Proc. R. Sco. Lond. A, 447,
155- 174 (1994).

50[9] Z. A. Munir and U. Anselmi-Tamburini, “Self-propagating Exothermic Reactions: theSynthesis of High-temperature Materials by Combustion”, North-Holland (1989).
[10] J. A. Pojman, V. M. Ilyashenko and A. M. Khan, “Spin mode instabilities in propagating fronts of polymerization”, Physica D84, 260-268 (1995).

[11] G. Sivashinsky, “On Spinning Propagation of Combustion Waves”, SIAM J. Appl.
Math. 40 (3), 432-438 (1981).
[12] T. Yamamoto, “Chemical kinetics in self-propagating high-temperature
Master’s thesis, Ryukoku University, (1995) (in Japanese).

$\mathrm{s}\mathrm{y}\mathrm{n}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{e}\mathrm{s}^{)}’$,