Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis

1. Introduction

The purpose of this book chapter is to analyze from a numerical point of view a reaction diffusion mathematical modelling of titanium carbide combustion synthesis from a mixture of titanium and carbide reactive powders thanks to self-propagating high temperature synthesis process. This modelling results in the coupling between a nonlinear parabolic equation expressing the enthalpy balance of the system with radiative boundary conditions and a nonlinear differential equation describing the exothermic chemical reaction in the system. An another multiphysics coupling was analyzed in [3]. This Self-propagating High temperature Synthesis (SHS) process was discovered in 1965 by Merzhanov [7], [8] and uses the energy released by the exothermic reaction itself to ensure its self-propagation inside the material after a localized heat supply has been performed for several seconds on the surface of the solid mixture. The stoichiometric solid mixture is made of several kinds of reactive powders. We analyze in this book chapter, the influence of radiative boundary conditions related to the heat supply over the ignition and eventual propagation of a combustion front inside the material.

Four sections are used to present our numerical simulation work. Section two presents the governing equations of the modelling. Section three outline the main aspects of the numerical scheme. Section four analyzes and discusses the main numerical simulation results of the combustion synthesis process. A conclusion summarizes the results that were obtained.

2. Mathematical modelling

This section describes the main features of the modelling which expresses the coupling between a reaction-diffusion written for the enthalpy balance and a differential equation written for the exothermic chemical kinetics. SHS (Self-propagation High-temperature Synthesis) process is a condensed phase process which converts a mixture of powders into an end product. In this paper we consider the synthesis of titanium carbide in solid phase ( $\left(s\right)$ subscript) thanks to the exothermic reaction $\mathrm{T}{\mathrm{i}}_{\left(s\right)}+{\mathrm{C}}_{\left(s\right)}\to \mathrm{Ti}{\mathrm{C}}_{\left(s\right)}$. We assume that we have a stoichiometric and isotropic mixture of compacted powders of titanium and carbide. Let us denote by $T=T\left(M,t\right)$ ( resp. $Ϝ=Ϝ\left(M,t\right)$) the temperature ( resp. conversion rate) at point $M\nabla \mathrm{\chi }$ at time $t$. The system is composed of the fraction $Ϝ$ of titanium carbide TiC end product and the fraction $(1-Ϝ)$ of remaining powders of titanium and carbide.

A first order exothermic kinetics is used to describe the synthesis of titanium carbide, hence a single variable, the conversion rate $Ϝ\nabla \left[\mathrm{0,1}\right]$ of the reaction is defined. When $Ϝ=0$ ( resp. $Ϝ=1$) the reaction has not started ( resp. is ended). The velocity of the reaction may have or not a cutoff temperature $T_s=1166$K, which corresponds to the first allotropic phase transition of titanium $Ϗ$, Ti${}_{Ϗ}$ to titanium $ϐ$, Ti${}_{ϐ}$. An Arrhenius type equation is used

$$\frac{dϜ}{\mathit{dt}}=k_d\left(T\right)\left(1-Ϝ\right),Ϝ(\mathrm{.},0)=0\mathrm{.}$$E1

The velocity constant $k_d\left(T\right)$ follows an Arrhenius law in temperature such that

$$k_d\left(T\right)=k_d^{\mathrm{*}}{T}^d\exp \left(-\frac{E^{\mathrm{*}}}{R\mathrm{.}T}\right)\mathrm{.}$$E2

In both cases, (with or without cut-off temperature), $k_d^{\mathrm{*}}$ is a frequency factor in $s^{-1}$, $E^{\mathrm{*}}$ is the activation energy of the reaction in $J/\mathit{mol}$, $R=8.314J/\left(\mathit{mol}\mathrm{.}K\right)$ is the perfect gas constant and $d\nabla \left[\mathrm{0,2}\right]$ is the degree, which is not necessarily an integer. Such an expression is commonly used in literature. It is worth mentioning that the activation temperature $T^{\mathrm{*}}=E^{\mathrm{*}}/R=4000\mathrm{}K$ is high and gives an idea of the stiffness of the reaction-diffusion/kinetics coupling.

2.2. Heat transfer modelling

This subsection presents successively the expressions used for the heat capacity and the thermal conductivity. It then analyzes the enthalpy balance that expresses the coupling between the heat transfer by thermal diffusion with the exothermic kinetics and provides the boundary conditions to close the mathematical modelling of the system. We also define and compute the adiabatic temperature of the exothermic kinetics. Due to the high temperatures involved in the process, up to 3000K, two phase changes occur

• $T=T_{Ϗ,ϐ}=1166$$\mathrm{T}{\mathrm{i}}_{Ϗ}$$\mathrm{T}{\mathrm{i}}_{ϐ}$$L_{Ϗ,ϐ}$17$f_{\mathit{Ϗϐ}}\nabla \left[\mathrm{0,1}\right]$$ϐ$
• $T=T_{\mathit{sl}}=1943$$L_{\mathit{sl}}$17$f_{\mathit{sl}}\nabla \left[\mathrm{0,1}\right]$

2.2.1. Equation used for the heat capacity

The mass heat capacity at constant pressure is a function of temperature, conversion rate and porosity. Assuming that the mass of the system remains constant during the process, we can therefore apply a linear mixing law between the heat capacity of reactants $C_{p_{\mathit{Ti}}}\left(T\right),C_{p_C}\left(T\right)$ and the heat capacity of the product $C_{p_{\mathit{TiC}}}\left(T\right)$ weighted by the conversion rate of the reaction $Ϝ$ to obtain

$$C_p\left(T,Ϝ\right)=Ϝ\mathrm{.}{C}_{p_{{}_{\mathit{TiC}}}}\left(T\right)+\left(1-Ϝ\right)\mathrm{.}\frac{M_{\mathit{Ti}}\mathrm{.}{C}_{p_{{}_{\mathit{Ti}}}}\left(T\right)+M_C\mathrm{.}{C}_{p_C}\left(T\right)}{M_{\mathit{TiC}}},$$E3

where $M_{\mathit{Ti}}$,$M_C$ et $M_{\mathit{TiC}}$ denote respectively the molar mass of titanium, carbide and titanium-carbide. The expression given by Eq.(3) is also used in [5],[6].

2.2.2. Equation used for the thermal conductivity $ϙ\left(T,Ϝ,f_{\mathit{sl}}\right)$

Thermal conductivity $ϙ\left(T,Ϝ,f_{\mathit{sl}}\right)$ in $J\mathrm{.}{m}^{-1}\mathrm{.}{K}^{-1}\mathrm{.}{s}^{-1}$ is a key parameter governing the propagation of the combustion front. Following [14] we define the effective thermal conductivity $ϙ\left(T,Ϝ,f_{\mathit{sl}}\right)$ by

$$ϙ\left(T,Ϝ,f_{\mathit{sl}}\right)={ϙ}_{\mathit{cond}}\left(T,Ϝ,f_{\mathit{sl}}\right)+{ϙ}_{\mathit{fus}}\left(f_{\mathit{sl}}\right)+{ϙ}_{\mathit{rad}}\left(T\right)+{ϙ}_{\mathit{conv}}\left(T\right),$$E4

where

• ${ϙ}_{\mathit{cond}}\left(T,Ϝ\right)$${ϙ}_{\mathrm{Ti}+\mathrm{C}}$${ϙ}_{\mathit{TiC}}\left(T_a\right)$

$${ϙ}_{\mathit{cond}}\left(T,Ϝ\right)=\left\{{ϙ}_{\mathit{Ti}+C}\left(T\right)+\left({ϙ}_{\mathit{TiC}}\left(T\right)-{ϙ}_{\mathit{Ti}+C}\left(T\right)\right)\mathrm{.}\sqrt{Ϝ}\right\}\mathrm{.}f\left(T\right),$$E5

while the expression of function $f\left(T\right)$ is given in [17]. Expression used for ${ϙ}_{\mathit{TiC}}\left(T\right),{ϙ}_{\mathit{Ti}+C}\left(T\right)$ are given in [13].

• ${ϙ}_{\mathit{fus}}\left(f_{\mathit{sl}}\right)=f_{\mathit{sl}}\left(1-Ϝ\right)$$T=T_{\mathit{sl}}=1943$$f_{\mathit{sl}}$$\left(1-Ϝ\right)$
• ${ϙ}_{\mathit{rad}}\left(T\right)$$d_p$$ϓ$$p$

$${ϙ}_{\mathit{rad}}\left(T\right)=4\mathrm{}\mathit{ϓϡ}{T}^3{d}_p{p}^{-3}\mathrm{.}$$E6

• ${ϙ}_{\mathit{conv}}\left(T\right)$${ϙ}_{\mathit{air}}\left(T\right)$

$${ϙ}_{\mathit{conv}}\left(T\right)={ϙ}_{\mathit{air}}\left(T\right)p^{-3}\mathrm{.}$$E7

2.2.3. Enthalpy balance

The enthalpy balance expressing the local conservation of internal energy can be written as

$$ϟ\frac{\mathit{dh}}{\mathit{dt}}=∆\mathrm{.}\left(ϙ\mathrm{}∆T\right),$$E8

where $ϟ$ in $\mathit{kg}\mathrm{.}{m}^{-3}$ is the density of the system, $h$ the mass enthalpy and $ϙ$ the thermal conductivity. Temperature $T$ and conversion rate $Ϝ$ are two thermodynamical variables of the problem. Moreover two scalar variables $f_{\mathit{Ϗϐ}}$ and $f_{\mathit{sl}}$, representing the fraction of solid and liquid phases are also used for the description of the mass enthalpy of the system. We have therefore a set of four independent variables $T,Ϝ,f_{\mathit{Ϗϐ}},f_{\mathit{sl}}$ that describe the evolution of the system. The mass enthalpy is defined by

$$h\left(T,Ϝ,f_{\mathit{Ϗϐ}},f_{\mathit{sl}}\right)=Ϝ{\left({\mathrm{\beta }}_{r}H\right)}_{T_a}+{\cup }_{T_a}^T\mathrm{‎}{C}_p\left(\varpi ,Ϝ\right)d\varpi +f_{\mathit{Ϗϐ}}{L}_{\mathit{Ϗϐ}}+f_{\mathit{sl}}{L}_{\mathit{sl}}\mathrm{.}$$E9

The computation of the partial derivatives $\frac{\complement h}{\complement T}$, $\frac{\complement h}{\complement Ϝ}$, $\frac{\complement h}{\complement f_{\mathit{Ϗϐ}}}$ and $\frac{\complement h}{\complement f_{\mathit{sl}}}$ leads to

$$\frac{\complement h}{\complement T}=C_p\left(T,Ϝ\right),\mathrm{}\frac{\complement h}{\complement f_{\mathit{Ϗϐ}}}=L_{\mathit{Ϗϐ}},\mathrm{}\frac{\complement h}{\complement f_{\mathit{sl}}}=L_{\mathit{sl}},\frac{\complement h}{\complement Ϝ}={\left({\mathrm{\beta }}_{r}H\right)}_{T_a}+{\cup }_{T_a}^T\mathrm{‎}\frac{\complement }{\complement Ϝ}\left[C_p\left(\varpi ,Ϝ\right)\right]\varpi \mathrm{.}$$E10

Thanks to the definition of $C_p\left(T,Ϝ\right)$ given by Eq.(3), we obtain

$$\frac{\complement h}{\complement Ϝ}={\left({\mathrm{\beta }}_{r}H\right)}_{T_a}+{\cup }_{T_a}^T\mathrm{‎}\left[C_{p_{{}_{\mathit{TiC}}}}\left(\varpi \right)-\frac{M_{\mathit{Ti}}\mathrm{.}{C}_{p_{{}_{\mathit{Ti}}}}\left(\varpi \right)+M_C\mathrm{.}{C}_{p_C}\left(\varpi \right)}{M_{\mathit{TiC}}}\right]d\varpi ,$$E11

where the integral $I\left(T\right)$ is computed with the trapezoidal rule.

$$I\left(T\right)={\cup }_{T_a}^T\mathrm{‎}\left[C_{p_{\mathit{TiC}}}\left(\varpi \right)-\frac{{M_{\mathit{TiC}}}_{p_{\mathit{Ti}}}\left(\varpi \right)+{M_C}_{p_C}\left(\varpi \right)}{M_{\mathit{TiC}}}\right]d\varpi \mathrm{.}$$E12

Fig.1 represents the evolution of $\frac{I\left(T\right)}{{\mathrm{\beta }}_{r}H}$for $T\nabla \left[\mathrm{300,3500}\right]$. It is observed that ${\max }_{300≣T≣3500}\left|\frac{I\left(T\right)}{{\left({\mathrm{\beta }}_{r}H\right)}_{T_a}}\right|≣2.510^{-5}$. The contribution of ${\cup }_{T_a}^T\mathrm{‎}\frac{\complement }{\complement Ϝ}\left[C_p\left(\varpi ,Ϝ\right)\right]d\varpi$ to $\frac{\complement h}{\complement Ϝ}$ is therefore neglected.

Finally the enthalpy balance for the system is written by

$$ϟ\mathrm{.}\frac{\complement h_s\left(T,Ϝ\right)}{\complement t}=∆\mathrm{.}\left(ϙ\mathrm{}∆T\right)+ϟ\mathrm{.}{\left(-{\mathrm{\beta }}_{r}H\right)}_{T_a}\mathrm{.}\frac{\complement Ϝ}{\complement t}-ϟ\mathrm{.}{L}_{\mathit{Ϗϐ}}\mathrm{.}\frac{\complement f_{\mathit{Ϗϐ}}}{\complement t}-ϟ\mathrm{.}{L}_{\mathit{sl}}\mathrm{.}\frac{\complement f_{\mathit{sl}}}{\complement t},$$E13

where

$$h_s\left(T,Ϝ\right)={\cup }_{T_a}^T\mathrm{‎}{C}_p\left(\varpi ,Ϝ\right)\varpi \mathrm{.}$$E14

Radiative boundary conditions are defined over $\complement \mathrm{\chi }$ for the system and take the form

$$-ϙ\frac{\complement T}{\complement n}=ϓ\mathrm{.}ϡ\mathrm{.}\left(T^4-T_{{\mathrm{\propto }}_{\complement \mathrm{\chi }}}^4\right),$$E15

where $ϓ=0.7$ is the emissivity of the material while $ϡ$ is the Stefan-Boltzman constant. The value of $T_{{\mathrm{\propto }}_{\complement \mathrm{\chi }}}$ depends on the boundary $\complement \mathrm{\chi }$.

The initial condition given on $\mathrm{\chi }$ states that the sample is at room temperature.

$$T\left(M,0\right)=T_a\mathrm{.}$$E16

Figure 1.

Evolution of$\frac{I\left(T\right)}{{\mathrm{\beta }}_{r}H}$for T ∈[300, 3500] K.

2.2.4. Adiabatic temperature of the reaction $T_{\mathit{ad}}$

The adiabatic temperature $T_{\mathit{ad}}$ of the reaction, or flame temperature, is the maximum temperature when the system is adiabatic. Titanium-carbide is obtained from titanium and carbide elements at temperature $T_a$. The enthalpy of the reaction corresponds to the heat of formation of TiC at temperature $T_a$, hence ${\left({\mathrm{\beta }}_{r}H\right)}_{T_a}={\mathrm{\beta }}_f{H}^0{\left(\mathit{TiC}\right)}_{T_a}$. Considering a thermodynamical path composed of the two steps:

1. Titanium and carbide reactants are transformed into titanium-carbide, -TiC- at $T_a$,

2. End product TiC is heated up from $T_a$ to $T_{\mathit{ad}}$.

One obtains that

$${\mathrm{\beta }}_f{H}^0{\left(\mathit{TiC}\right)}_{T_a}+{\cup }_{T_a}^{T_{\mathit{ad}}}\mathrm{‎}{C}_p\left(\mathrm{TiC}\right)\mathit{dT}=0\mathrm{.}$$E17

We use Janaf tables [11] for the expression of the heat capacity as a function of temperature. When $T_a=298K$ (room temperature value), we solve numerically the integral equation Eq.(17), and obtain that the value of the adiabatic temperature is $T_{\mathit{ad}}=2900K$. According to [17], we express $-{\mathrm{\beta }}_{r}H=\frac{\mathrm{\beta }{H}_f}{M_{\mathit{TiC}}}$.

2.3. Mathematical modelling of SHS process

The purpose of the numerical simulation is to determine at each point $M\nabla \mathrm{\chi }$, and at each time $t$, temperature field $T(M,t)$, conversion rate $Ϝ(M,t)$, solid fraction $f_{\mathit{Ϗϐ}}(M,t)$ and liquid fraction $f_{\mathit{sl}}(M,t)$ that will permit to determine spatial and temporal location of the reaction ignition. The mathematical modelling given by Eq.(1),Eq.(2) and Eq.(13),Eq.(14),Eq.(15) is well posed and will be solved numerically by the methods described in the next section.

3. Numerical discretization scheme

In this section we present the main features of the fully implicit finite-volume discretization scheme such as discrete maximum principle for both the reaction-diffusion and the differential equation discretization [1]. Error estimates are given in [2]. The iterative solution of penta-diagonal sparse matrix for 2D computations and hepta-diagonal sparse matrix for 3D computations is done thanks to SSOR and SIP methods. A fixed point technique is used to solve the coupled nonlinear system. A first order linearization of the exponential Arrhenius term is used which enhances the diagonal dominance of the matrix and accelerates the iterative SSOR/SIP solvers. We have used the enthalpy method to compute the phase change for which a detailed description can be found in [10]. This will be omitted from now, since we will mainly focuss our attention on the construction of the numerical scheme for the reaction-diffusion equation on a non uniform structured mesh and discard formally the specifics of the enthalpy method.

3.1. One-dimensional finite-volume discretization

We present the set of discrete equations that arise from the implicit finite-volume discretization of the enthalpy balance (reaction-diffusion equation) and the exothermic chemical reaction (differential equation).

3.1.1. Numerical discretization of the enthalpy balance

We integrate the enthalpy balance over a space-time finite volume ${\mathrm{\chi }}_i\times \left[t_n,t_{n+1}\right]$, where ${\mathrm{\chi }}_i=[r_i,r_{i+1}]$ is a control volume. A cell-centered approximation is used. The unknown $T_i^n$ denotes the mean value of $T(x,t)$ at time $t_n$ over ${\mathrm{\chi }}_i$. $m_i={\cup }_{r_i}^{r_{i+1}}\mathrm{‎}{r}^g\mathit{dr}$ is the discrete "mass" of control volume ${\mathrm{\chi }}_i$.

$${\cup }_{t_n}^{t_{n+1}}\mathrm{‎}{\cup }_{{\mathrm{\chi }}_i}\mathrm{‎}ϟ\mathrm{.}\frac{\complement h_s\left(T,Ϝ\right)}{\complement t}d\mathrm{\chi }\mathit{dt}={\cup }_{t_n}^{t_{n+1}}\mathrm{‎}{\cup }_{{\mathrm{\chi }}_i}\mathrm{‎}\left(∆\left(ϙ\left(T,Ϝ,f_{\mathit{sl}}\right)∆\right)+\frac{ϟ\mathrm{.}\mathrm{\beta }{H}_f}{M_{\mathit{TiC}}}k\left(T\right)\left(1-Ϝ\right)\right)d\mathrm{\chi }\mathit{dt}\mathrm{.}$$E18

Backward Euler implicit scheme is used for the temporal discretization since it is unconditionally stable, robust and well adapted for the discretization of such problems. $\mathrm{\beta }t$ is the time-step.

We give only for internal control volumes, the finite-volume discretization of the enthalpy balance in a one-dimensional cartesian ($g=0$), cylindrical ($g=1$) and spherical ($g=2$) coordinate system

$$m_iϟ\frac{h_s\left(T_i^{n+1},{Ϝ}_i^{n+1}\right)-h_s\left(T_i^n,{Ϝ}_i^n\right)}{\mathrm{\beta }t}={ϙ}_{i+\frac{1}{2}}^{n+1}{r}_{i+1}^g\frac{T_{i+1}^{n+1}-T_i^{n+1}}{d_{i+\frac{1}{2}}}-{ϙ}_{i-\frac{1}{2}}^{n+1}{r}_i^g\frac{T_i^{n+1}-T_{i-1}^{n+1}}{d_{i-\frac{1}{2}}}+m_i\frac{ϟ\mathrm{\beta }{H}_f}{M_{\mathit{TiC}}}{k}_d\left(T^{n+1}\right)\left(1-{Ϝ}_i^{n+1}\right)\mathrm{.}$$E19

Here $d_{i+\frac{1}{2}}$ is the distance between the center of two adjacent control volumes, and denoting by $h_i=r_{i+1}-r_i$ the length of control volume ${\mathrm{\chi }}_i$, then $d_{i+\frac{1}{2}}=(h_{i+1}+h_i)/2$.

As pointed out by [9], in order to ensure the conservativity of the heat flux at the interface between two adjacents control volumes, ${ϙ}_{i+\frac{1}{2}}^{n+1}$ is computed thanks to the following harmonic mean formula:

$${ϙ}_{i+\frac{1}{2}}^{n+1}=\frac{h_i+h_{i+1}}{\frac{h_{i+1}}{{ϙ}_{i+1}^{n+1}}+\frac{h_i}{{ϙ}_i^{n+1}}},$$E20

where ${ϙ}_i^{n+1}=ϙ\left(T_i^{n+1},{Ϝ}_i^{n+1},{f_{\mathit{sl}}}_i^{n+1}\right)$. A decoupled iterative solution of the nonlinear system is achieved thanks to the fixed point method. A first-order linearization of both the stiff Arrhenius term, and the enthalpy term is done. This reinforces the strictly dominant three-diagonal matrix that is inverted at each step of the non-linear solver thanks to the direct $(Ϗ,ϐ)$ Gauss algorithm, i.e. the TDMA algorithm [9]. The numerical solution cost of the three-diagonal linear system by this method increases linearly with the number of unknowns.

3.1.2. Numerical discretization of the exothermic kinetics

We integrate the differential equation over a space-time finite volume ${\mathrm{\chi }}_i\times \left[t_n,t_{n+1}\right]$. A cell-centered approximation is used.

$${\cup }_{t_n}^{t_{n+1}}\mathrm{‎}{\cup }_{{\mathrm{\chi }}_i}\mathrm{‎}\frac{dϜ}{\mathit{dt}}={\cup }_{t_n}^{t_{n+1}}\mathrm{‎}{\cup }_{{\mathrm{\chi }}_i}\mathrm{‎}{k}_d\left(T\right)\left(1-Ϝ\right)\mathrm{.}$$E21

The backward Euler scheme, first order accurate fully implicit scheme is used. Denoting by ${Ϝ}_i^n$ the mean value at time $t_n$ over ${\mathrm{\chi }}_i$ of $Ϝ(x,t)$, we obtain the following discrete equation

$${Ϝ}_i^{n+1}=\frac{{Ϝ}_i^n+\mathrm{\beta }t{k}_d\left(T_i^{n+1}\right)}{1+\mathrm{\beta }t{k}_d\left(T_i^{n+1}\right)}\mathrm{.}$$E22

It was proved in [1] that the numerical scheme is $L^{\mathrm{\propto }}$ stable, moreover for each time index $n$ and for all control volume index $i=1,‧,I\mathrm{}:\mathrm{}{Ϝ}_i^n\nabla \left[\mathrm{0,1}\right]$. For a given control volume index $i$, the time sequence $({Ϝ}_i^n{)}_n$ is increasing, hence the discrete finite-volume approximation mimics the behavior of the continuous solution [2].

3.2. Two-dimensional finite-volume discretization

We will not detail the set of discrete equations for the chemical kinetics, since it is a straightforward extension from the 1D case. We now give the discrete equations related to the finite-volume approximation of the enthalpy balance over a space-time finite volume ${\mathrm{\chi }}_{i,j}\times \left[t_n,t_{n+1}\right]$. Each rectangular structured non uniform control volume ${\mathrm{\chi }}_{i,j}=[r_i,r_{i+1}]\times [z_j,z_{j+1}]$, has surface $m_{i,j}$ defined by :

$$m_{i,j}={\cup }_{{\mathrm{\chi }}_{i,j}}\mathrm{‎}d\mathrm{\chi }={\cup }_{z_j}^{z_{j+1}}\mathrm{‎}\left({\cup }_{r_i}^{r_{i+1}}\mathrm{‎}{r}^g\right)\mathit{dz}\mathrm{}=\frac{{\left(r_{i+1}\right)}^{g+1}-{\left(r_i\right)}^{g+1}}{g+1}\mathrm{.}\left(z_{j+1}-z_j\right)\mathrm{.}$$E23

The enthalpy balance is written for 2D cartesian geometry $(g=0)$ or 2D cylindrical geometry $(g=1)$.

$\frac{m_{i,j}}{\mathrm{\beta }t}\mathrm{.}ϟ\mathrm{.}\left(h_s\left(T_{i,j}^{n+1},{Ϝ}_{i,j}^{n+1}\right)-h_s\left(T_{i,j}^n,{Ϝ}_{i,j}^n\right)\right)={ϙ}_{i+\frac{1}{2},j}^{n+1}{r}_{i+1}^g\frac{T_{i+1,j}^{n+1}-T_{i,j}^{n+1}}{d_{i+\frac{1}{2}}^r}-{ϙ}_{i-\frac{1}{2},j}^{n+1}{r}_i^g\frac{T_{i,j}^{n+1}-T_{i-1,j}^{n+1}}{d_{i-\frac{1}{2}}^r}+{ϙ}_{i,j+\frac{1}{2}}^{n+1}\frac{T_{i,j+1}^{n+1}-T_{i,j}^{n+1}}{d_{j+\frac{1}{2}}^z}-{ϙ}_{i,j-\frac{1}{2}}^{n+1}\frac{T_{i,j}^{n+1}-T_{i,j-1}^{n+1}}{d_{j-\frac{1}{2}}^z}+m_{i,j}\mathrm{.}\frac{ϟ\mathrm{.}\mathrm{\beta }{H}_f}{M_{\mathit{TiC}}}k\left(T_{i,j}^{n+1}\right)\left(1-{Ϝ}_{i,j}^{n+1}\right)\mathrm{.}$

The expression used for ${ϙ}_{i+\frac{1}{2},j}^{n+1}$, ${ϙ}_{i-\frac{1}{2},j}^{n+1}$, ${ϙ}_{i,j+\frac{1}{2},j}^{n+1}$, ${ϙ}_{i,j-\frac{1}{2}}^{n+1}$ is a straightforward adaptation of Eq.(20). We use the same decoupled iterative solution of the nonlinear system as in the one-dimensional case. A first-order linearization of both the stiff Arrhenius term, and the enthalpy term is done. This reinforces the strictly dominant sparse penta-diagonal matrix that is inverted at each step of the non-linear solver. It is known that the more strictly dominant a matrix is, the faster iterative solver such as the SSOR, (successive over relaxation method) converges.

3.3. Three-dimensional finite-volume discretization

We now give the discrete equations related to the finite-volume approximation of the enthalpy balance over a space-time finite volume ${\mathrm{\chi }}_{i,j,k}\times \left[t_n,t_{n+1}\right]$. Each parallelepipedic structured non uniform control volume ${\mathrm{\chi }}_{i,j,k}=[x_i,x_{i+1}]\times [x_j,y_{j+1}]\times [z_k,z_{k+1}]$, has volume $m_{i,j,k}={\cup }_{{\mathrm{\chi }}_{i,j,k}}\mathrm{‎}d\mathrm{\chi }$.

$$\frac{m_{i,j,k}}{\mathrm{\beta }t}\mathrm{.}ϟ\mathrm{.}\left(h_s\left(T_{i,j,k}^{n+1},{Ϝ}_{i,j,k}^{n+1}\right)-h_s\left(T_{i,j,k}^n,{Ϝ}_{i,j,k}^n\right)\right)={ϙ}_{i+\frac{1}{2},j,k}^{n+1}\frac{T_{i+1,j,k}^{n+1}-T_{i,j,k}^{n+1}}{d_{i+\frac{1}{2}}^x}-{ϙ}_{i-\frac{1}{2},j,k}^{n+1}\frac{T_{i,j,k}^{n+1}-T_{i-1,j,k}^{n+1}}{d_{i-\frac{1}{2}}^x}+{ϙ}_{i,j+\frac{1}{2},k}^{n+1}\frac{T_{i,j+1,k}^{n+1}-T_{i,j,k}^{n+1}}{d_{j+\frac{1}{2}}^y}-{ϙ}_{i,j-\frac{1}{2},k}^{n+1}\frac{T_{i,j,k}^{n+1}-T_{i,j-1,k}^{n+1}}{d_{j-\frac{1}{2}}^y}+{ϙ}_{i,j,k+\frac{1}{2}}^{n+1}\frac{T_{i,j,k+1}^{n+1}-T_{i,j,k}^{n+1}}{d_{k+\frac{1}{2}}^z}-{ϙ}_{i,j,k-\frac{1}{2}}^{n+1}\frac{T_{i,j,k}^{n+1}-T_{i,j,k-1}^{n+1}}{d_{k-\frac{1}{2}}^z}+m_{i,j,k}\mathrm{.}\frac{ϟ\mathrm{.}\mathrm{\beta }{H}_f}{M_{\mathit{TiC}}}k\left(T_{i,j,k}^{n+1}\right)\left(1-{Ϝ}_{i,j,k}^{n+1}\right)\mathrm{.}$$E24

The expression used for ${ϙ}_{i+\frac{1}{2},j,k}^{n+1}$, ${ϙ}_{i-\frac{1}{2},j,k}^{n+1}$, ${ϙ}_{i,j+\frac{1}{2},j,k}^{n+1}$, ${ϙ}_{i,j-\frac{1}{2},k}^{n+1}$, ${ϙ}_{i,j,k+\frac{1}{2}}^{n+1}$, ${ϙ}_{i,j,k-\frac{1}{2}}^{n+1}$ is a straightforward adaptation of Eq.(20). We use the same numerical procedure as the one used for the two-dimensional case. It is worth mentioning that a sparse strictly dominant hepta-diagonal matrix is inverted at each step of the non-linear solver.

3.4. Numerical implementation

Due to space constraints, we only describe key aspects of our numerical software entitled Hephaïstos; a toolbox library for multidimensional numerical computation of combustion synthesis problems written in C/C++. Detailed documented results related to the implementation on various single core, multi-core and many-core architectures are given in [4]. Auto-parallelization and openmp based speedup results on core i7 quad-core using gnu gcc, sun studio and open64 compilers are reported. A similar study is done on cuda-core using nvidia nvcc compiler [4]. Efforts were done to achieve the highest possible performance on processors that use memory cache hierarchy such as MIPS R1X00 processors and INTEL Core i7 processors. As an example, all loops invariant quantities are pre-computed and stored into one-dimensional arrays. Dynamically allocated multidimensional arrays are not used, because of pointer aliasing problems. We prefer to convert such arrays into one-dimensional arrays. A storage similar to CRS(Compact Row Storage) is used that takes into account the penta (hepta)-diagonal sparsity of the matrices that are iteratively inverted by SSOR and SIP methods. The postprocessing of the software’s results saved to vtk format is done thanks to mayavi software for temporal snapshot and Paraview for interactive analysis and movie production.

4. 1D/2D/3D numerical study

In this section we analyze the ignition and propagation of the combustion front during titanium carbide-TiC- combustion synthesis. It is worth mentioning that the propagation of the combustion wave is either stable i.e. at constant velocity or oscillatory and can be determined according to [16] by the computation of the Zeldovich number $\mathit{Ze}$ defined by $\mathit{Ze}=\frac{\left(T_{\mathit{ad}}-T_a\right)E^{\mathrm{*}}}{R_{\mathit{ad}}^2}$. $T_{\mathit{ad}}\left[K\right]$ is the adiabatic temperature of the reaction, $T_a\left[K\right]$ is the ambiant temperature and $\left(k_d^{\mathrm{*}}\left[s^{-1}\right],E^{\mathrm{*}}[\mathit{kJ}/\mathit{mol}]\right)$ characterizes the first order exothermic kinetics. $R$ is the gaz constant. Experimental observations were reported in [8] to sustain this theoretical analysis. Moreover in the constant velocity case resp.( oscillatory case i.e. periodic oscillation of the velocity of the solid flame propagation around a mean velocity), the synthesised titanium-carbide product has uniform resp.( non-uniform) physical properties. For $(k_0,E^{\mathrm{*}})=(3800\mathrm{}{s}^{-1},30.0\mathrm{}\mathit{kcal}\mathrm{.}\mathit{mo}{l}^{-1})$, the Zeldovich number $\mathit{Ze}=4.52<Z{e}_c=2\left(2+\sqrt{5}\right)$, therefore the solid flame propagation is stable.

4.1. 1D numerical study

This subsection presents unpublished results related to the induction time ${Ϣ}_{\mathit{ind}}$ in planar/cylindrical/spherical geometry such as one depicted in Fig.(21).

The boundary conditions used to analyze the combustion front propagation in the reactive mixture of length $R_e$ are

$$-ϙ\frac{\complement T}{\complement n}\left(0,t\right)=\mathit{ϓϡ}\left(T{\left(0,t\right)}^4-{\left(T_{f_{x=0}}\right)}^4\right),-ϙ\frac{\complement T}{\complement n}\left(R_e,t\right)=\mathit{ϓϡ}\left(T{\left(R_e,t\right)}^4-{\left(T_{f_{x=R_e}}\right)}^4\right)\mathrm{.}$$E25

4.1.1. Tools used to analyse the numerical simulation results

We define the induction time, starting point and ending time, three notions that will be heavily used in this section.

• The induction time ${Ϣ}_{\mathit{ind}}$ is the required time for the build-up of a steady-state propagation regime of the chemical reaction inside the material. It is defined as the first instant for which there exists $x\nabla \mathrm{\chi }$ such that $Ϝ(x,{Ϣ}_{\mathit{ind}})>0$. It depends on various parameters such as ignition temperature, heat capacity, thermal conductivity, pre-heating time, mass density.

• Starting of the reaction is $x_{\mathit{ind}}>0$ such that $Ϝ(x_{\mathit{ind}},{Ϣ}_{\mathit{ind}})=0.5$,

• Ending time of the reaction is time ${Ϣ}_{\mathit{end}}>0$ such that $⇿x\nabla \mathrm{\chi },Ϝ(x,{Ϣ}_{\mathit{end}})=1$,

• Thickness of reaction zone. Assuming that the reaction starts at ${Ϣ}_{\mathit{ind}}$, then we can follow the evolution of the combustion front propagation in the material thanks to the computation of spatial and temporal evolution of points $x_M^{0.01}\left(t\right)$, $x_M^{0.50}\left(t\right)$ et $x_M^{0.99}\left(t\right)$ such that $Ϝ\left(x_M^{0.01}\right(t),t)=0.01$, $Ϝ\left(x_M^{0.50}\right(t),t)=0.50$ et $Ϝ\left(x_M^{0.99}\right(t),t)=0.99$,

• A each instant $t$, $x_{\mathit{max}}\left(t\right)$ denotes the location which has the maximum temperature in domain $\mathrm{\chi }$,

• The synthesis temperature characterizes the thermal history of the material from the point of view of the kinetics and is defined for $x\nabla \mathrm{\chi }$ by

$$T_{\mathit{syn}}\left(x\right)={\cup }_0^{+}\mathrm{‎}T(x,t)\frac{dϜ}{\mathit{dt}}(x,t)\mathrm{}\mathit{dx}\mathrm{.}$$E26

It was defined and used in [2].

A meaningful numerical simulation is presented in Fig.(2) and in Fig.(3). We observe that the "thickness" of the reaction-zone defined by $x_M^{0.99}\left(t\right)-x_M^{0.01}\left(t\right)$ is nearly constant. Moreover, the slope of the three curves is nearly the same in the steady state regime. This confirms that a stable propagation anticipated thanks to the computation of the Zeldovich number is effectively obtained. Nevertheless, the temporal evolution of $x_{\mathit{max}}\left(t\right)$ gives an information about the thermal history. Fig.(2) and Fig.(3) represent the temporal evolution of $x_{\mathit{max}}\left(t\right)$, $x(Ϝ=0.01,t)\mathrm{}x(Ϝ=0.50,t)$ and $x(Ϝ=0.99,t)$ for wall temperature $T_{f_{r=0}}=1600$ K. Phase changes are taken ( resp. not taken ) into account in Fig.(3) (resp in Fig.(2)). One observes on both Fig.(2) and Fig.(3) that $x_{\mathit{max}}\left(t\right)$ is equal to zero during preheating time, moreover taking into account phase changes as seen on Fig.(3) modifies the velocity of $x_{\mathit{max}}\left(t\right)$. When the reaction starts, the front propagates, and $x_{\mathit{max}}\left(t\right)$ evolves differently than $x_M^{0.01}\left(t\right),x_M^{0.99}\left(t\right)$. Moreover when the front reaches the extremity of the sample, then $x_{\mathit{max}}\left(t\right)$ increases rapidly and soon after decreases because the synthesis is finished and the radiative heat losses contribute significantly to the decreasing of the temperature field inside the material.

Figure 2.

Temporal evolution of $x_{\mathit{max}}\left(t\right)$ and $x(Ϝ=0.01,t)\mathrm{}x(Ϝ=0.50,t)\mathrm{}x(Ϝ=0.99,t)$ for wall temperature $T_{f_{r=0}}=1600$ K. Phase changes are not taken into account.

Figure 3.

Temporal evolution of $x_{\mathit{max}}\left(t\right)$ and $x(Ϝ=0.01,t)\mathrm{}x(Ϝ=0.50,t)\mathrm{}x(Ϝ=0.99,t)$ for wall temperature $T_{f_{r=0}}=1600$ K. Phase changes are taken into account.

4.1.2. Combustion front propagation

Fig.(4) shows that the propagation of the stiff combustion front is stable, since each profile in the steady state propagation regime are equally distant from each other.

Figure 4.

Temporal evolution of $T(x,t)$ at equally spaced instants. Phase change is not taken into account.

Fig.(5) shows that the energy released by the exothermic kinetics is nearly constant during the propagation of the combustion wave inside the material.

Figure 5.

Temporal evolution of $\frac{ϟ\mathrm{\beta }{H}_f}{M_{\mathit{TiC}}}{k}_d\left(T\right)\left(1-{Ϝ}_i\right)$ at equally spaced instants. Phase change is not taken into account.

4.1.3. Contribution of furnace temperature

A sample of length $R_e=1$cm is placed inside a furnace maintained at temperature $T_f\nabla \left[\mathrm{800,2400}\right]$K. Fig.(6) shows that both induction time ${Ϣ}_{\mathit{ind}}$ and ending time ${Ϣ}_{\mathit{end}}$ increase exponentially when $T_f$ decreases, whenever phase change is taken into account or not.

Figure 6.

Evolution of ${Ϣ}_{\mathit{ind}}$ et ${Ϣ}_{\mathit{end}}$ as a function of temperature furnace $T_{f_{r=0}}$.

Fig.(7) shows that the evolution of $x_{\mathit{ind}}\left(t\right)$ is nearly similar whatever the furnace temperature $T_f$ is. In fact, each curve can be translated from each other, so the behaviour is nearly independent from $T_f$, except from small values for which a slight curvature effect is observed. A similar conclusion can be drawn upon analyzing Fig.(8) which represents the time evolution of $T_{\mathit{max}}\left(t\right)$. It is also worth mentioning that a sharp peak is observed on each curve. An explanation comes from the fact that when the combustion front, at high temperature ends its propagation, it reaches the cold extremity, so an intense heat transfer occurs.

Figure 7.

Temporal evolution of $x_{\mathit{ind}}\left(t\right)$ position, for various values of the furnace temperature $T_{f_{r=0}}\nabla \left[\mathrm{800,1200,1600,2000,2400}\right]K$.

Figure 8.

Time evolution of $T_{\mathit{max}}\left(t\right)$ for several values of temperature furnace $T_{f_{r=0}}\nabla \left[\mathrm{800,1200,1600,2000,2400}\right]$.

Fig.(9) shows the same behaviour for $x_{\mathit{max}\left(t\right)}$ whatever $T_f$ is. The following analysis is proposed.

• A pre-heating stage, for which thanks to the heat supply at $x=0$ (radiative boundary condition), there is a gradual increase in temperature. But $x_{\mathit{max}}\left(t\right)=0$, so the temperature is below the ignition temperature.

• A combustion front propagation step $x_{\mathit{max}}\left(t\right)$ increases linearly with respect to time, so a constant velocity propagation of TiC synthesis is observed.

• An acceleration of the propagation which comes from the influence of the radaitive boundary condition at $x=R_e$. During a small time interval, the hot spot reaches the ending extremity of the sample.

• A cooling stage by thermal conduction, since the combustion synthesis is finished. The hot spot moves rapidly towards the center of the material thanks to the cooling. It is also observed that the velocity of the hot spot is a function of $T_f$.

Finally in the steady state propagation regime, each curve $x_{\mathit{max}}\left(t\right)$ can be translated from each other. The local heat supply induced by the boundary condition at $x=0$ doesn’t modify the characteristics of the combustion front propagation such as it’s velocity. The time evolution of $T_{\mathit{max}\left(t\right)}$ depicted in Fig.(8) can be analyzed similarly.

Figure 9.

Time evolution of $x_{\mathit{max}}\left(t\right)$ for different values of furnace temperature $T_{f_{r=0}}$.

Spatial distribution for different values of $T_{f_{r=0}}$ for temperature ( resp. conversion rate) at $t=1$s, is represented on Fig.(10). Using a suitable scaling, it is observed that their shape and stiffness is similar.

Figure 10.

Spatial distribution of temperature profile at $t=1$s, for each wall temperature $T_f\nabla \left\{800,\mathrm{}1200,\mathrm{}1600,\mathrm{}2000,\mathrm{}2400\right\}\mathrm{K}$.

Figure 11.

Spatial distribution of conversion rate profile at $t=1$s, for each wall temperature $T_f\nabla \left\{800,\mathrm{}1200,\mathrm{}1600,\mathrm{}2000,\mathrm{}2400\right\}\mathrm{K}$.

4.1.4. Energy released/absorbed by the boundary conditions

We analyze the time evolution of flux ${\mathrm{\tau }}_0\left(t\right)$, exchanged at $r=0$ between the exterior and the material, defined by ${\mathrm{\tau }}_0\left(t\right)=-\left(ϙ\frac{\complement T}{\complement n}\right)\left(0,t\right)=\mathit{ϓϡ}\left(T(0,t{)}^4-T_f^4\right)$. Three consecutive steps can be observed thanks to the analysis of ${\mathrm{\tau }}_0\left(t\right)$ represented by Fig.(12) and correlated to the time-evolution of $T(0,t)$

Figure 12.

Time evolution of the energy released/removed by the radiative boundary conditions when phase change is taken into account or not.

1. $T(0,t)$$T_a$$T_f$${\mathrm{\tau }}_0\left(t\right)$
2. $T(0,t)$$T_f$$T_{\mathit{ad}}$$T_f<T_{\mathit{ad}}$${\mathrm{\tau }}_0\left(t\right)$
3. $T(0,t)$$T_f$${\mathrm{\tau }}_0\left(t\right)$

A similar analysis can be done for ${\mathrm{\tau }}_{R_e}\left(t\right)=-\left(ϙ\frac{\complement T}{\complement n}\right)\left(R_e,t\right)=\mathit{ϓϡ}\left(T(R_e,t{)}^4-T_f^4\right)$. It is worth mentioning that due to the thermal shock observed when the "hot" solid combustion front reaches the "cold" boundary, the amplitude ${\mathrm{\tau }}_{R_e}\left(t\right)$ is an order of magnitude higher than ${\mathrm{\tau }}_0\left(t\right)$. Taking into account phase change doesn’t modify the observation done thanks to Fig.(12).

4.1.5. Contribution of the cut-off temperature to the ignition and propagation

We assume that the chemical kinetics has a cut-off temperature $T_s$ equal to the first phase transition $T{i}_{Ϗ}\to T{i}_{ϐ}$ temperature, $T_{\mathit{Ϗϐ}}=1166$K. We observe on Fig.(13), that the discrepancy ${Ϣ}_{\mathit{ign}}^{T_s=0}-{Ϣ}_{\mathit{ign}}^{T_s=T_{\mathit{Ϗϐ}}}$ increases significantly when $T_f$ decreases whenever the phase change is taken into account or not. In practice, ${\left({Ϣ}_{\mathit{ign}}^{T_s=0}\right)}_{T_f=900\mathrm{}K}={\left({Ϣ}_{\mathit{ign}}^{T_s=T_{\mathit{Ϗϐ}}}\right)}_{T_f=1600\mathrm{}K}$. A similar conclusion is drawn on Fig.(14) for the difference ${Ϣ}_{\mathit{end}}^{T_s=0}-{Ϣ}_{\mathit{end}}^{T_s=T_{\mathit{Ϗϐ}}}$.

Below the cut-off temperature, the kinetics is not active, therefore the modelling reduces to a non-linear diffusion equation. The main phenomena is the pre-heating of the material. When the temperature reaches locally the ignition temperature for a certain amount of time, the chemical reaction starts.

Figure 13.

Evolution of ${Ϣ}_{\mathit{ind}}$ as a function of $t_{f_{r=0}}$ for a kinetics with/without a cut-off temperature. Phase change taken into account ($\mathit{PC}=1$), or not ($\mathit{PC}=0$).

Figure 14.

Evolution of ${Ϣ}_{\mathit{ind}}$ as a function of $t_{f_{r=0}}$ for a kinetics with/without a cut-off temperature. Phase change taken into account ($\mathit{PC}=1$), or not ($\mathit{PC}=0$).

We analyze the contribution of the cut-off temperature to the spatial stiffness of the combustion front through the temperature profile on Fig.(15) and conversion rate profile on Fig.(16). In both cases, the spatial resolution of the numerical discretisation scheme is satisfactory, and we observe that temperature and conversion rate profiles are stiffer when a cut-off temperature is taken into account.

Figure 15.

Spatial temperature profile at $t=2$s, for chemical kinetic without cutoff temperature $T_s=0$, and with cutoff temperature $T_s=T_{\mathit{Ϗϐ}}$.

Figure 16.

Spatial conversion rate profile at $t=2$s, for chemical kinetic without cutoff temperature $T_s=0$, and with cutoff temperature $T_s=T_{\mathit{Ϗϐ}}$.

4.1.6. Analysis of a double front propagation

We consider the situation when we heat identically both extremities of the cylinder with $T_{f_{r=0}}=T_{f_{r=R_e}}=T_f=2400\mathrm{}$K. Obviously spatial profiles are similar and symmetrical with respect to $x=R_e/2$ as can be seen on Fig.(17)-(18).

Figure 17.

Spatial temperature profiles at several consecutive instants.

Figure 18.

Spatial conversion rate profiles at several consecutive instants.

Moreover, when the two fronts meet each other at $x=R_e/2$, a liquid phase may appear. Temperature at the center of the material is significantly greater than the adiabatic temperature of the reaction represented in Fig.(19). The stiffness of the conversion rate in the double front case versus the single front case is also noticed in Fig.(20).

Figure 19.

Temporal evolution of $T(R_e,t)$ for single and double fronts propagation.

Figure 20.

Temporal evolution of $Ϝ(R_e,t)$ for single and double fronts propagation.

4.1.7. Contribution of the geometry to the induction time

Fig.(21) shows the evolution of $\log \left({Ϣ}_{\mathit{ind}}\right)\left(T_f\right)$ for slab, cylindrical and spherical geometry. The three profiles are similar, nevertheless it is observed that ${Ϣ}_{\mathit{ind}}^{\mathrm{slab}}>{Ϣ}_{\mathit{ind}}^{\mathrm{cyl}}>{Ϣ}_{\mathit{ind}}^{\mathrm{sph}}$, for each value of $T_f$. Curvature effects may explain this result.

Figure 21.

Evolution of the $\log \left({Ϣ}_{\mathit{ind}}\right)\left(T_f\right)$ for slab geometry $g=0$, cylindrical geometry $g=1$ and spherical geometry $g=2$.

4.1.8. Contribution of the kinetics to the induction time

We analyze the contribution of the exponent $d$ defined in Eq.(2) to the induction time ${Ϣ}_{\mathit{ind}}$ and the ending time ${Ϣ}_{\mathit{end}}$. In order to obtain a better precision for ${Ϣ}_{\mathit{ind}}\left(d\right)$, we use fractional exponents $d\nabla \left[\mathrm{0,2}\right]$. Moreover to be able to compare the results obtained for various values of $k_d\left(T\right)$, we perform a normalization of the pre-exponential factor such that its value remain the same when $T=T_{\mathit{ad}}$ (the adiabatic temperature). We write therefore the equality between factors $k_0\left(T\right)=k_0^{\mathrm{*}}\mathit{exp}\left(-\frac{E^{\mathrm{*}}}{\mathit{RT}}\right)$ and $k_d\left(T\right)=k_d^{\mathrm{*}}{T}^d\mathit{exp}\left(-\frac{E^{\mathrm{*}}}{\mathit{RT}}\right)$ at $T=T_{\mathit{ad}}$

$$k_0^{\mathrm{*}}=k_d^{\mathrm{*}}\mathrm{.}{T}_{\mathit{ad}}^d\mathrm{.}$$E27

Fig.(22) represents the temperature dependance of $k_d\left(T\right)$. Ten numerical simulations are done for five equally-spaced values of $d\nabla \left[\mathrm{0,2}\right]$. For each value of $d$ we take into account or not the phase change and use $T_{f_{r=0}}=1600$K.

Figure 22.

Temperature dependance of $k_d\left(T\right)$ for various values of exponent $d\nabla \left[\mathrm{0,2}\right]$.

We point out that

• Each curve $k_d\left(T\right)$ is an increasing function of temperature $T$ and has the same concavity,

• When $T≣T_{\mathit{ad}}$, then $k_d(T{)}_{d>0}$ is below $k_0\left(T\right)$,

• When $T>T_{\mathit{ad}}$, then $k_d\left(T\right)>>k_0\left(T\right)$. This imply that the heat released by the exothermic kinetics is significantly increasing when $d$ increases. Temperature "overshoot" of $T_{\mathit{max}}\left(t\right)$ can be observed on Fig.(30).

• We conclude that the velocity of the combustion flame is higher when $d=0$ than when $d>0$, and the temperature obtained is super-adiabatic as seen on Fig.(30).

Fig.(23) shows that the increasing rate of the induction time ${Ϣ}_{\mathit{ind}}$ is super-linear, while it is linear for the ending time ${Ϣ}_{\mathit{end}}$ when degree $d$ increases. It appears independent from the eventual phase change contribution.

Figure 23.

Comparison of induction time ${Ϣ}_{\mathit{ind}}$ and ending time ${Ϣ}_{\mathit{end}}$, when phase change is taken into account (PC=1) or not (PC=0), for several values of exponent $d\nabla \left[\mathrm{0,2}\right]$, with $T_{f_{r=0}}=1600$ K and $T_{f_{r=R_e}}=300$ K.

It is worth mentioning that all spatial temperature profiles represented on Fig.(24) have a discontinuity of the time-derivative when $T=T_{\mathit{sl}}$. This is explained by the contribution ${ϙ}_{\mathit{fus}}\left(f_{\mathit{sl}}\right)$ of the titanium melting to the thermal conductivity given by Eq.(4).

Figure 24.

Spatial distribution of temperature field $T(x,t=5\mathrm{}s)$ for different values of exponent $d$, when phase change is taken into account (PC=1) or not (PC=0) with $T_{f_{r=0}}=1600$ K and $T_{f_{r=R_e}}=300$ K.

Conversion rate profiles represented on Fig.(25) are sharper when phase change is taken into account.

Figure 25.

Spatial distribution of conversion rate field $Ϝ(x,t=5\mathrm{}s)$ for different values of exponent $d$, when phase change is taken into account (PC=1) or not (PC=0) with $T_{f_{r=0}}=1600$ K and $T_{f_{r=R_e}}=300$ K.

The spatial distribution of the synthesis temperature $T_{\mathit{syn}}(x,t=30\mathrm{}s)$ always increases when $d$ increases, whether phase change is taken into account or not as seen on Fig.(26).

Figure 26.

Spatial distribution of synthesis temperature $T_{\mathit{syn}}(x,t=30\mathrm{}s)$ for different values of exponent $d$, when phase change is taken into account (PC=1) or not (PC=0) with $T_{f_{r=0}}=1600$ K and $T_{f_{r=R_e}}=300$ K.

Fig.(27) shows that the time evolution of $x_{\mathit{max}}\left(t\right)$ has a similar shape for each value of $d$ phase change is taken into account or not. Considering high values of $d$ imply a significant delay for ${Ϣ}_{\mathit{ind}}$. Moreover the spatial distribution of the heat released by the kinetics, as seen on Fig.(28) for $t=5$s presents a maximum for $d=0$, and a minimum for $d=2$. This is correlated to the high values of induction time ${Ϣ}_{\mathit{ind}}$ when high values of $d$ are used. Phase change contributes significantly to the time evolution of the front’s position $x_{Ϝ=1/2}\left(t\right)$ as observed on Fig.(29). Without phase change, the evolution is nearly independent from the value of $d$ since the slope of $x_{Ϝ=1/2}\left(t\right)$ is nearly the same. Time evolution of the front’s maximum temperature $T_{\mathit{max}}\left(t\right)$ for different values of exponent $d$ as seen on Fig.(30) is in fact nearly independent from $d$ when a suitable time-translation is performed whenever phase change is taken into account or not.

Figure 27.

Time evolution of $x_{\mathit{max}}\left(t\right)$ for different values of exponent $d$ when phase change is taken into account (PC=1) or not (PC=0) with $T_{f_{r=0}}=1600$ K and $T_{f_{r=R_e}}=300$ K.

Figure 28.

Spatial distribution of energy $E(x,t)$ released by the exothermic kinetics for different values of exponent $d$ when phase change is taken into account (PC=1) or not (PC=0) with $T_{f_{r=0}}=1600$ K and $T_{f_{r=R_e}}=300$ K.

Figure 29.

Time evolution of the front’s position $x_{Ϝ=1/2}\left(t\right)$ for different values of exponent $d$ when phase change is taken into account (PC=1) or not (PC=0) with $T_{f_{r=0}}=1600$ K and $T_{f_{r=R_e}}=300$ K.

Figure 30.

Time evolution of the front’s maximum temperature $T_{\mathit{max}}\left(t\right)$ for different values of exponent $d$ when phase change is taken into account (PC=1) or not (PC=0) with $T_{f_{r=0}}=1600$ K and $T_{f_{r=R_e}}=300$ K.

4.2. 2D numerical study

This subsection computes the propagation of a dual radial/longitudinal front and a single longitudinal front for various values of the temperature furnace in which the cylindrical sample is placed. The main interest of this modelling with respect to the 1D case is to analyze the contribution of the lateral heat losses over the ignition and propagation of the combustion front inside the cylindrical sample. The boundary conditions are defined over $\complement \mathrm{\chi }={\mathrm{\alpha }}_{r=0}\cap {\mathrm{\alpha }}_{r=R_e}\cap {\mathrm{\alpha }}_{z=0}\cap {\mathrm{\alpha }}_{z=H_e}$ by

$$\left\{\begin{array}{l}⇿\left(r,z\right)\nabla {\mathrm{\alpha }}_{r=0}:\left(-ϙ\frac{\complement T\left(0,z,t\right)}{\complement n_{r=0}^e}\right)=0,\\ ⇿\left(r,z\right)\nabla {\mathrm{\alpha }}_{r=R_e}:\left(-ϙ\frac{\complement T\left(R_e,z,t\right)}{\complement n_{r=R_e}^e}\right)=ϓ\mathrm{.}ϡ\mathrm{.}\left(T{\left(R_e,z,t\right)}^4-{\left(T_{f_{z,r=R_e}}\left(t\right)\right)}^4\right),\\ ⇿\left(r,z\right)\nabla {\mathrm{\alpha }}_{z=0}:\left(-ϙ\frac{\complement T\left(r,0,t\right)}{\complement n_{z=0}^e}\right)=ϓ\mathrm{.}ϡ\mathrm{.}\left(T{\left(r,0,t\right)}^4-{\left(T_{f_{r,z=0}}\left(t\right)\right)}^4\right),\\ ⇿\left(r,z\right)\nabla {\mathrm{\alpha }}_{z=H_e}:\left(-ϙ\frac{\complement T\left(r,H_e,t\right)}{\complement n_{z=H_e}^e}\right)=ϓ\mathrm{.}ϡ\mathrm{.}\left(T{\left(r,H_e,t\right)}^4-{\left(T_{f_{r,z=H_e}}\left(t\right)\right)}^4\right)\mathrm{.}\end{array}\right.$$E28

The choice of the temperature triple $\left(T_{f_{z,r=R_e}}\left(t\right),\mathrm{}{T}_{f_{r,z=0}}\left(t\right),\mathrm{}{T}_{f_{r,z=H_e}}\left(t\right)\right)$ characterizes the way the heat supplied on the exterior surface of the cylinder (radius $R_e=2$cm, height $H_e=2$cm) will induce the propagation of the combustion front inside the cylinder, initially at room temperature $T=T_a=300\mathrm{}K$. Each computation will analyze the phenomena during $t=10\mathrm{}s$. The computational grid is uniform along $z$-axis and iso-volume along $r$-axis.

4.2.1. Reaction-diffusion vs thermal diffusion

We choose $\left(T_{f_{z,r=R_e}}\left(t\right),\mathrm{}{T}_{f_{r,z=0}}\left(t\right),\mathrm{}{T}_{f_{r,z=H_e}}\left(t\right)\right)=\left(300\mathrm{}K,1600\mathrm{}K,300\mathrm{}K\right)$ for both simulations. The thermal diffusion case means that the kinetics is not active, i.e. $k_d=0$, so the initial mixture of reactive powders is not transformed. At each point $(r,z)$ of the sample and at each instant $t>0$, we observe on Fig.(31) that $T_a≣T(r,z,t)≣T_f$, when the exothermic kinetics is not taken into account, moreover the numerical solution fulfills a discrete maximum principle. A similar principle is observed when the kinetics is taken into account, but the time evolution is significantly different because of the sharp rise in temperature when the synthesis starts.

Figure 31.

Time evolution of $T(\mathrm{0,0},t)$ et $T(R_e,0,t)$.

Moreover $T(R_e,0,t)<T(\mathrm{0,0},t)$ because of the radiative heat losses at $r=R_e$.

4.2.2. Single, dual longitudinal, dual longitudinal/single radial front

These three cases were considered and analyzed in detail in [2].

4.2.3. Dual radial-longitudinal front

We consider the case where $\left(T_{f_{z,r=R_e}}\left(t\right),\mathrm{}{T}_{f_{r,z=0}}\left(t\right),\mathrm{}{T}_{f_{r,z=H_e}}\left(t\right)\right)=\left(2400\mathrm{}K,2400\mathrm{}K,300\mathrm{}K\right)$. At the same time a longitudinal front is moving from bottom to top, while a radial front moves from the exterior surface of the cylinder to the center, as seen on the temperature distribution at time $t=3$s represented by Fig.(32). The conversion rate $Ϝ(r,z,3)$ has a similar shape. The energy released by the exothermic process is nearly uniformly distributed on $\mathrm{\chi }$ as seen on Fig.(33). A similar conclusion can be drawn upon the shape of the synthesis temperature $T_{\mathit{syn}}(r,z,5)$ represented on Fig.(34), except along the line where the two fronts are joining themselves.

Figure 32.

Temperature distribution $T(r,z,3)$ for $(r,z)\nabla \mathrm{\chi }$.

Figure 33.

Integral of the energy released by the kinetics during the process at each $(r,z)\nabla \mathrm{\chi }$.

According to [2], in order to analyze the kind of propagation, we assume given a temperature $\mathrm{\zeta }\nabla \left[\mathrm{300,3000}\right]$ and a simulation time $t_f>0$. We define for each point $x\nabla \mathrm{\chi }$, the thermal history time $t_{\mathrm{\zeta }}\left(x\right)$ which corresponds to the total time for which a given point $x\nabla \mathrm{\chi }$ remains at a temperature greater or equal to $\mathrm{\zeta }$ and gives an indication of the thermal history of the material as a function of the coupling between the exothermic reaction and the thermal diffusion. More precisely, $t_{\mathrm{\zeta }}\left(x\right)$ determines if the energy involved in the reaction-diffusion process is uniformly distributed or not in $\mathrm{\chi }$ and what is the average temperature of the process. Fig.(35) represents such time distribution at $t=5$s for $\mathrm{\zeta }=1800$K. Fig.(36) represents such time distribution at $t=5$s for $\mathrm{\zeta }=2400$K. As in the previous figures, it is influenced by the propagation of the combustion wave. It is worth mentioning that on these two figures the time distribution is nearly equal, from up to 4.16s (resp. 4.70) for 1800K (resp. 2400K), and that they can be superposed.

Figure 34.

Synthesis temperature distribution $T_{\mathit{syn}}(r,z,5)$ for $(r,z)\nabla \mathrm{\chi }$.

Figure 35.

Time distribution of the thermal history time $T_{1800K}(r,z,5)$ for $(r,z)\nabla \mathrm{\chi }$.

Figure 36.

Time distribution of the thermal history time $T_{2400K}(r,z,5)$ for $(r,z)\nabla \mathrm{\chi }$.

4.3. 3D numerical study

This subsection accounts for the heterogeneity of the radiative boundary conditions to analyze the ignition and propagation of the combustion front in a cubic sample. We assume that the heat supplied to the six faces of the cubic sample is different, in the sense that the wall temperature that contributes to the preheating of each face of the cubic sample is different. This will induce a specific transient spatial pattern of the combustion front. A meaningful snapshot of the transient temperature profile is depicted in Fig.(37). It is clearly observed that the shape of the front’s propagation inside the cube is influenced by the boundary conditions. If a regular propagation is required, providing heat supply over one single face of the cube is enough to obtain quickly titanium-carbide. So combining a different heat supply system for each face of the cube appears complicated from an experimental point of view. Numerical simulation show that the pattern of the propagation in this case is more complicated than in the single heat supply case and requires more computational power to get finely resolved. An analysis, not presented here due space constraints, and using the same methodology as in the previous 2D case show that the same conclusions can be drawn for $T_{2400K}(x,y,z,\mathrm{.}),\mathrm{}{T}_{\mathit{syn}}(x,y,z,\mathrm{.})\mathrm{.}$

Figure 37.

Temperature distribution $T(x,y,z,3)$ for $(x,y,z)\nabla \mathrm{\chi }$.

5. Conclusions and perspectives

In this book chapter, we have presented a multidimensional modelling for the numerical computation of ignition and propagation of combustion fronts during combustion synthesis of ceramic materials such as titanium-carbide. A detailed computational study was done in 1D slab/cylindrical/spherical geometry, 2D cylindrical/cartesian and 3D cartesian geometry to analyze the influence of the radiative boundary conditions over the induction time. The radiation contribution to the thermal conductivity was taken into account and the sensitivity of the induction time to several parameters such as the kinetics, wall temperature, phase change was carefully analyzed. Our numerical software Hephaïstos was presented. It implements an implicit finite-volume scheme for which error estimates, discrete maximum principles were reported and used to ensure the consistency of the numerical results. This modelling study was done for titanium carbide TiC. It can be applied to other ceramic materials such as silicium carbide-SiC-. Moreover in order to analyze the combustion front propagation for high Zeldovich number cases, an adaptive finite-volume scheme is required. A new monotonicity preserving refinement/derefinement conservative algorithm has been designed for multidimensional computations in various coordinate systems. This algorithm maintains the structured topology of the mesh. It is currently under implementation.