**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

©2012 Aoufi and Damamme, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

( (*s*) subscript) thanks to the exothermic reaction Ti(s) + C(s) → TiC(s). We assume that we have a stoichiometric and isotropic mixture of compacted powders of titanium and carbide. Let us denote by *T* = *T* (*M*, *t*) (*resp. ξ* = *ξ* (*M*, *t*)) the temperature (*resp. conversion rate*) at point *M* ∈ Ω 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.

### **2.1. Exothermic kinetics modelling**

A first order exothermic kinetics is used to describe the synthesis of titanium carbide, hence a single variable, the conversion rate *ξ* ∈ [0, 1] 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 *Ts* = 1166K, which corresponds to the first allotropic phase transition of titanium *α*, Ti*<sup>α</sup>* to titanium *β*, Ti*β*. An Arrhenius type equation is used

$$\begin{aligned} \frac{d\mathfrak{F}}{dt} &= k\_d \begin{pmatrix} T \end{pmatrix} \begin{pmatrix} 1 - \mathfrak{F} \end{pmatrix}, \\ \mathfrak{F}(.,0) &= 0. \end{aligned} \tag{1}$$

and the heat capacity of the product *Cp*TiC (*T*) weighted by the conversion rate of the reaction

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 503

where *M*Ti,*M*<sup>C</sup> et *M*TiC denote respectively the molar mass of titanium, carbide and

Thermal conductivity *λ* (*T*, *ξ*, *fsl*) in *J*.*m*−1.*K*−1.*s*−<sup>1</sup> is a key parameter governing the propagation of the combustion front. Following [14] we define the effective thermal

• *λcond* (*T*, *ξ*) is the thermal conductivity of the system which is a non-linear weighting between the thermal conductivity of the reactants *λ*Ti+<sup>C</sup> and the product *λTiC*(*Ta*)

*λ*Ti+C(*T*) + (*λ*TiC(*T*) − *λ*Ti+C(*T*)).

while the expression of function *f*(*T*) is given in [17]. Expression used for

• *λf us* (*fsl*) = *fsl* (1 − *ξ*) accounts for the appearance of a liquid phase when temperature *T* = *Tsl* = 1943K, i.e. the temperature for the fusion of titanium. In this case, the fraction *fsl* will increase from 0 to 1. Moreover this temperature is reached thanks to the exothermic kinetics while a fraction (1 − *ξ*) of titanium has already been consumed during the kinetics.

• *λrad* (*T*) expresses the contribution of the radiative heat transfer between the grains and depends on the diameter *dp* of the titanium/carbide particules, their emissivity *�*, their porosity *p*. The expression used is taken from [18]. It was reported in [15] that this term

• *λconv* (*T*) expresses the contribution of gaz that are present in the porous system. *λair* (*T*)

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

*ρ dh*

This expression represents the fraction of fused titanium that has not yet reacted.

*λ* (*T*, *ξ*, *fsl*) = *λcond* (*T*, *ξ*, *fsl*) + *λf us* (*fsl*) + *λrad* (*T*) + *λconv* (*T*), (4)

*M*Ti.*Cp*Ti (*T*) + *M*C.*Cp*<sup>C</sup> (*T*) *M*TiC

> *ξ*

*<sup>λ</sup>rad* (*T*) = <sup>4</sup> *ε σ <sup>T</sup>*<sup>3</sup> *dp <sup>p</sup>*<sup>−</sup>3. (6)

*λconv* (*T*) = *λair* (*T*) *p*<sup>−</sup>3. (7)

*dt* <sup>=</sup> <sup>∇</sup>.(*<sup>λ</sup>* <sup>∇</sup>*T*), (8)

, (3)

. *f* (*T*), (5)

*Cp* (*T*, *ξ*) = *ξ*.*Cp*TiC (*T*) + (1 − *ξ*).

*2.2.2. Equation used for the thermal conductivity λ* (*T*, *ξ*, *fsl*)

titanium-carbide. The expression given by Eq.(3) is also used in [5],[6].

*ξ* to obtain

where

conductivity *λ* (*T*, *ξ*, *fsl*) by

is taken from [13].

*2.2.3. Enthalpy balance*

according to [17]. It is given by

*λ*TiC(*T*), *λ*Ti+C(*T*) are given in [13].

*λcond* (*T*, *ξ*) =

contributes significantly to the velocity of the combustion wave.

The velocity constant *kd*(*T*) follows an Arrhenius law in temperature such that

$$k\_d\left(T\right) = k\_d^\* T^d \exp\left(-\frac{E^\*}{R.T}\right). \tag{2}$$

In both cases, (with or without cut-off temperature), *k*∗ *<sup>d</sup>* is a frequency factor in *<sup>s</sup>*−1, *<sup>E</sup>*<sup>∗</sup> is the activation energy of the reaction in *J*/*mol*, *R* = 8.314*J*/(*mol*.*K*) is the perfect gas constant and *d* ∈ [0, 2] 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*∗ = *E*∗/*R* = 4000 *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


#### *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 *Cp*Ti (*T*), *Cp*<sup>C</sup> (*T*) and the heat capacity of the product *Cp*TiC (*T*) weighted by the conversion rate of the reaction *ξ* to obtain

$$\mathbb{C}\_{p}\left(T,\xi\right) = \xi.\mathbb{C}\_{p\_{\text{TC}}}\left(T\right) + \left(1 - \xi\right).\frac{M\_{\text{Ti}}.\mathbb{C}\_{p\_{\text{Ti}}}\left(T\right) + M\_{\text{C}}.\mathbb{C}\_{p\_{\text{C}}}\left(T\right)}{M\_{\text{TiC}}},\tag{3}$$

where *M*Ti,*M*<sup>C</sup> et *M*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 λ* (*T*, *ξ*, *fsl*)

Thermal conductivity *λ* (*T*, *ξ*, *fsl*) in *J*.*m*−1.*K*−1.*s*−<sup>1</sup> is a key parameter governing the propagation of the combustion front. Following [14] we define the effective thermal conductivity *λ* (*T*, *ξ*, *fsl*) by

$$
\lambda\left(T,\mathfrak{f}\_{\prime}f\_{\rm sl}\right) = \lambda\_{\rm cond}\left(T,\mathfrak{f}\_{\prime}f\_{\rm sl}\right) + \lambda\_{\rm fus}\left(f\_{\rm sl}\right) + \lambda\_{\rm rad}\left(T\right) + \lambda\_{\rm conv}\left(T\right), \tag{4}
$$

where

2 Numerical Simulations

( (*s*) subscript) thanks to the exothermic reaction Ti(s) + C(s) → TiC(s). We assume that we have a stoichiometric and isotropic mixture of compacted powders of titanium and carbide. Let us denote by *T* = *T* (*M*, *t*) (*resp. ξ* = *ξ* (*M*, *t*)) the temperature (*resp. conversion rate*) at point *M* ∈ Ω at time *t*. The system is composed of the fraction *ξ* of titanium carbide **TiC** end

A first order exothermic kinetics is used to describe the synthesis of titanium carbide, hence a single variable, the conversion rate *ξ* ∈ [0, 1] 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 *Ts* = 1166K, which corresponds to the first allotropic phase transition of

*<sup>d</sup> <sup>T</sup><sup>d</sup>* exp

activation energy of the reaction in *J*/*mol*, *R* = 8.314*J*/(*mol*.*K*) is the perfect gas constant and *d* ∈ [0, 2] 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*∗ = *E*∗/*R* = 4000 *K*

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

• at *T* = *Tα*,*<sup>β</sup>* = 1166K, the allotropic phase transition of Ti*<sup>α</sup>* to Ti*β*, characterized by its mass

• at *T* = *Tsl* = 1943K, titanium melting is characterized by its mass latent heat *Lsl*,[17], and

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 *Cp*Ti (*T*), *Cp*<sup>C</sup> (*T*)

is high and gives an idea of the stiffness of the reaction-diffusion/kinetics coupling.

 <sup>−</sup> *<sup>E</sup>*<sup>∗</sup> *R*.*T* 

*dt* <sup>=</sup> *kd* (*T*) (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*), (1)

. (2)

*<sup>d</sup>* is a frequency factor in *<sup>s</sup>*−1, *<sup>E</sup>*<sup>∗</sup> is the

product and the fraction (1 − *ξ*) of remaining powders of titanium and carbide.

titanium *α*, Ti*<sup>α</sup>* to titanium *β*, Ti*β*. An Arrhenius type equation is used

*dξ*

*ξ*(., 0) = 0.

The velocity constant *kd*(*T*) follows an Arrhenius law in temperature such that

*kd* (*T*) = *k*<sup>∗</sup>

In both cases, (with or without cut-off temperature), *k*∗

involved in the process , up to 3000K, two phase changes occur

*fsl* ∈ [0, 1] is the fraction of liquid phase.

*2.2.1. Equation used for the heat capacity*

latent heat *Lα*,*β*,[17], and *fαβ* ∈ [0, 1] is the fraction of *β* phase,

**2.2. Heat transfer modelling**

**2.1. Exothermic kinetics modelling**

• *λcond* (*T*, *ξ*) is the thermal conductivity of the system which is a non-linear weighting between the thermal conductivity of the reactants *λ*Ti+<sup>C</sup> and the product *λTiC*(*Ta*) according to [17]. It is given by

$$
\lambda\_{\text{cond}}\left(T,\xi\right) = \left\{\lambda\_{\text{Ti}+\mathbb{C}}\left(T\right) + \left(\lambda\_{\text{Ti}\mathbb{C}}\left(T\right) - \lambda\_{\text{Ti}+\mathbb{C}}\left(T\right)\right) . \sqrt{\xi}\right\} . f\left(T\right), \tag{5}
$$

while the expression of function *f*(*T*) is given in [17]. Expression used for *λ*TiC(*T*), *λ*Ti+C(*T*) are given in [13].


$$
\lambda\_{rad} \left( T \right) = 4 \,\varepsilon \,\sigma \, T^3 \, d\_p \, p^{-3}. \tag{6}
$$

• *λconv* (*T*) expresses the contribution of gaz that are present in the porous system. *λair* (*T*) is taken from [13].

$$
\lambda\_{conv}(T) = \lambda\_{air}(T) \ p^{-3}.\tag{7}
$$

#### *2.2.3. Enthalpy balance*

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

$$
\rho \frac{dh}{dt} = \nabla. \left(\lambda \,\nabla T\right),
\tag{8}
$$

#### 4 Numerical Simulations 504 Numerical Simulation – From Theory to Industry Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis <sup>5</sup>

where *ρ* in *kg*.*m*−<sup>3</sup> 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αβ* and *fsl*, 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αβ*, *fsl* that describe the evolution of the system. The mass enthalpy is defined by

$$\ln \left( T\_{\prime} \mathfrak{F}\_{\prime} f\_{\mathbf{a} \not\! \not\! \mathbf{f}} \, f\_{\mathrm{sl}} \right) = \mathfrak{F}\_{\prime} \left( \Delta\_{\prime} H \right)\_{T\_{\mathrm{d}}} + \int\_{T\_{\mathrm{d}}}^{T} \mathbb{C}\_{p} \left( \theta\_{\prime} \mathfrak{F}\_{\prime} \right) \, d\theta + f\_{\mathbf{a} \not\!\! \mathbf{f}} \, \mathcal{L}\_{\mathbf{a} \not\!\! \mathbf{f}} + f\_{\mathrm{sl}} \, \mathcal{L}\_{\mathrm{sl}}.\tag{9}$$

The computation of the partial derivatives *<sup>∂</sup><sup>h</sup> <sup>∂</sup><sup>T</sup>* , *∂h ∂ξ* , *<sup>∂</sup><sup>h</sup> ∂ fαβ* and *<sup>∂</sup><sup>h</sup> ∂ fsl* leads to

$$\frac{\partial \hbar}{\partial T} = \mathbb{C}\_p \left( T, \tilde{\xi} \right), \quad \frac{\partial \hbar}{\partial f\_{a\beta}} = \mathrm{L}\_{a\beta\nu} \quad \frac{\partial \hbar}{\partial f\_{s\ell}} = \mathrm{L}\_{s\ell\nu}$$

$$\frac{\partial \hbar}{\partial \tilde{\xi}} = (\Delta\_r H)\_{T\_\delta} + \int\_{T\_\delta}^T \frac{\partial}{\partial \tilde{\xi}} \left[ \mathrm{C}\_p \left( \theta, \tilde{\xi} \right) \right] \, d\theta. \tag{10}$$

300 1300 2300 3300 **T [K]**

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 505

The adiabatic temperature *Tad* 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 *Ta*. The enthalpy of the reaction corresponds to the heat

. Considering a

*Cp* (TiC) *dT* = 0. (17)

−3.0e−05

*2.2.4. Adiabatic temperature of the reaction Tad*

thermodynamical path composed of the two steps:

(ii) End product TiC is heated up from *Ta* to *Tad*.

<sup>Δ</sup>*rH* for *<sup>T</sup>* <sup>∈</sup> [300, 3500] K.

of formation of TiC at temperature *Ta*, hence (Δ*rH*)*Ta* <sup>=</sup> <sup>Δ</sup>*<sup>f</sup> <sup>H</sup>*<sup>0</sup> (TiC)*Ta*

<sup>Δ</sup>*<sup>f</sup> <sup>H</sup>*<sup>0</sup> (*TiC*)*Ta* <sup>+</sup>

*MTiC* .

**2.3. Mathematical modelling of SHS process**

(i) Titanium and carbide reactants are transformed into titanium-carbide, -TiC- at *Ta*,

 *Tad Ta*

We use Janaf tables [11] for the expression of the heat capacity as a function of temperature. When *Ta* = 298*K* (room temperature value), we solve numerically the integral equation Eq.(17), and obtain that the value of the adiabatic temperature is *Tad* = 2900*K*. According

The purpose of the numerical simulation is to determine at each point *M* ∈ Ω, and at each time *t*, temperature field *T*(*M*, *t*), conversion rate *ξ*(*M*, *t*), solid fraction *fαβ*(*M*, *t*) and liquid fraction *fsl*(*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.

**Figure 1.** Evolution of *<sup>I</sup>*(*T*)

One obtains that

to [17], we express <sup>−</sup>Δ*rH* <sup>=</sup> <sup>Δ</sup>*Hf*

−2.0e−05

−1.0e−05

0.0

1.0e−05

Thanks to the definition of *Cp* (*T*, *ξ*) given by Eq.(3), we obtain

$$\frac{\partial \hbar}{\partial \xi} = (\Delta\_{\rm r} H)\_{T\_d} + \int\_{T\_d}^{T} \left[ \mathbf{C}\_{p\_{\rm TC}} \left( \theta \right) - \frac{M\_{\rm Ti} \mathbf{C}\_{p\_{\rm Ti}} \left( \theta \right) + M\_{\rm C} \mathbf{C}\_{p\_{\rm C}} \left( \theta \right)}{M\_{\rm TiC}} \right] d\theta,\tag{11}$$

where the integral *I*(*T*) is computed with the trapezoidal rule.

$$I(T) = \int\_{T\_d}^{T} \left[ \mathbb{C}\_{p\_{\rm TIC}} \left( \theta \right) - \frac{M\_{\rm TIC} \mathbb{C}\_{p\_{\rm TIC}} \left( \theta \right) + M\_{\rm C} \mathbb{C}\_{p\_{\rm C}} \left( \theta \right)}{M\_{\rm TIC}} \right] d\theta. \tag{12}$$

Fig.1 represents the evolution of *<sup>I</sup>*(*T*) <sup>Δ</sup>*rH* for *<sup>T</sup>* <sup>∈</sup> [300, 3500]. It is observed that max 300≤*T*≤3500 *I* (*T*) (Δ*rH*)*Ta* <sup>≤</sup> 2.5 10<sup>−</sup>5. The contribution of *<sup>T</sup> Ta ∂ ∂ξ Cp* (*θ*, *ξ*) *dθ* to *∂h ∂ξ* is therefore neglected.

Finally the enthalpy balance for the system is written by

$$\rho.\frac{\partial \mathfrak{h}\_{\sf s}\left(T,\xi\right)}{\partial t}=\nabla.\left(\lambda\,\nabla T\right)+\rho.\left(-\Delta\_{\sf r}H\right)\_{T\_{\sf s}}.\frac{\partial \mathfrak{f}}{\partial t}-\rho.L\_{\sf a\sf \beta}.\frac{\partial f\_{\sf a\sf \beta}}{\partial t}-\rho.L\_{\sf s}.\frac{\partial f\_{\sf s\sf \sf}}{\partial t},\tag{13}$$

where

$$h\_{\mathfrak{s}}\left(T,\mathfrak{f}\right) = \int\_{T\_{\mathfrak{s}}}^{T} \mathbb{C}\_{p}\left(\theta,\mathfrak{f}\right) \,d\theta. \tag{14}$$

Radiative boundary conditions are defined over *∂*Ω for the system and take the form

$$-\lambda \, \frac{\partial T}{\partial n} = \varepsilon \, \sigma \, \left(T^4 - T\_{\infty \text{in}}^4\right),\tag{15}$$

where *�* = 0.7 is the emissivity of the material while *σ* is the Stefan-Boltzman constant. The value of *T*∞*∂*<sup>Ω</sup> depends on the boundary *∂*Ω.

The initial condition given on Ω states that the sample is at room temperature.

$$T\left(M\_\prime 0\right) = T\_\mathfrak{a}.\tag{16}$$

**Figure 1.** Evolution of *<sup>I</sup>*(*T*) <sup>Δ</sup>*rH* for *<sup>T</sup>* <sup>∈</sup> [300, 3500] K.

### *2.2.4. Adiabatic temperature of the reaction Tad*

The adiabatic temperature *Tad* 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 *Ta*. The enthalpy of the reaction corresponds to the heat of formation of TiC at temperature *Ta*, hence (Δ*rH*)*Ta* <sup>=</sup> <sup>Δ</sup>*<sup>f</sup> <sup>H</sup>*<sup>0</sup> (TiC)*Ta* . Considering a thermodynamical path composed of the two steps:

(i) Titanium and carbide reactants are transformed into titanium-carbide, -TiC- at *Ta*,

(ii) End product TiC is heated up from *Ta* to *Tad*.

One obtains that

4 Numerical Simulations

where *ρ* in *kg*.*m*−<sup>3</sup> 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αβ* and *fsl*, 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αβ*, *fsl* that describe the evolution of the

> *T Ta*

*<sup>∂</sup><sup>T</sup>* , *∂h ∂ξ* , *<sup>∂</sup><sup>h</sup> ∂ fαβ*

*∂ fαβ*

*∂ ∂ξ Cp* (*θ*, *ξ*)

 *T Ta*

*Cp* (*θ*, *ξ*) *dθ* + *fαβ Lαβ* + *fsl Lsl*. (9)

leads to

*dθ*. (10)

<sup>Δ</sup>*rH* for *<sup>T</sup>* <sup>∈</sup> [300, 3500]. It is observed that

*∂ fαβ*

*T* (*M*, 0) = *Ta*. (16)

*<sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>ρ</sup>*.*Lsl*.

*Cp* (*θ*, *ξ*) *dθ*. (14)

*dθ*, (11)

*dθ*. (12)

*∂ξ* is therefore

*<sup>∂</sup><sup>t</sup>* , (13)

*∂h*

*∂ fsl*

, (15)

and *<sup>∂</sup><sup>h</sup> ∂ fsl*

= *Lsl*,

*∂ fsl*

*M*TiC

<sup>=</sup> *<sup>L</sup>αβ*, *<sup>∂</sup><sup>h</sup>*

*Cp*TiC (*θ*) <sup>−</sup> *<sup>M</sup>*Ti.*Cp*Ti (*θ*) <sup>+</sup> *<sup>M</sup>*C.*Cp*<sup>C</sup> (*θ*)

*M*TiC

 *T Ta*

*∂ξ*

*∂ ∂ξ Cp* (*θ*, *ξ*) *dθ* to

*<sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>ρ</sup>*.*Lαβ*.

*Cp*TiC (*θ*) <sup>−</sup> *<sup>M</sup>*TiC *Cp*Ti (*θ*) <sup>+</sup> *<sup>M</sup>*<sup>C</sup> *Cp*<sup>C</sup> (*θ*)

system. The mass enthalpy is defined by

*T*, *ξ*, *fαβ*, *fsl*

The computation of the partial derivatives *<sup>∂</sup><sup>h</sup>*

*∂ξ* <sup>=</sup> (Δ*rH*)*Ta* <sup>+</sup>

*I* (*T*) =

Fig.1 represents the evolution of *<sup>I</sup>*(*T*)

  *∂h*

*∂h*

= *ξ* (Δ*rH*)*Ta* +

*<sup>∂</sup><sup>T</sup>* <sup>=</sup> *Cp* (*T*, *<sup>ξ</sup>*), *<sup>∂</sup><sup>h</sup>*

*∂ξ* <sup>=</sup> (Δ*rH*)*Ta* <sup>+</sup>

 *T Ta*

<sup>≤</sup> 2.5 10<sup>−</sup>5. The contribution of

*hs* (*T*, *ξ*) =

<sup>−</sup> *<sup>λ</sup> <sup>∂</sup><sup>T</sup>*

 *T Ta*

where *�* = 0.7 is the emissivity of the material while *σ* is the Stefan-Boltzman constant. The

*<sup>T</sup>*<sup>4</sup> <sup>−</sup> *<sup>T</sup>*<sup>4</sup> ∞*∂*<sup>Ω</sup> 

Radiative boundary conditions are defined over *∂*Ω for the system and take the form

*<sup>∂</sup><sup>n</sup>* <sup>=</sup> *<sup>ε</sup>*.*σ*.

The initial condition given on Ω states that the sample is at room temperature.

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>∇</sup>.(*<sup>λ</sup>* <sup>∇</sup>*T*) <sup>+</sup> *<sup>ρ</sup>*.(−Δ*rH*)*Ta* .

Thanks to the definition of *Cp* (*T*, *ξ*) given by Eq.(3), we obtain

where the integral *I*(*T*) is computed with the trapezoidal rule.

 *T Ta* 

Finally the enthalpy balance for the system is written by

*h* 

*∂h*

max 300≤*T*≤3500

neglected.

where

 

*ρ*.

*I* (*T*) (Δ*rH*)*Ta*

*∂hs* (*T*, *ξ*)

value of *T*∞*∂*<sup>Ω</sup> depends on the boundary *∂*Ω.

$$
\Delta\_f H^0 \left( \text{TiC} \right)\_{T\_d} + \int\_{T\_d}^{T\_{ad}} \mathbb{C}\_p \left( \text{TiC} \right) \, dT = 0. \tag{17}
$$

We use Janaf tables [11] for the expression of the heat capacity as a function of temperature. When *Ta* = 298*K* (room temperature value), we solve numerically the integral equation Eq.(17), and obtain that the value of the adiabatic temperature is *Tad* = 2900*K*. According to [17], we express <sup>−</sup>Δ*rH* <sup>=</sup> <sup>Δ</sup>*Hf MTiC* .

#### **2.3. Mathematical modelling of SHS process**

The purpose of the numerical simulation is to determine at each point *M* ∈ Ω, and at each time *t*, temperature field *T*(*M*, *t*), conversion rate *ξ*(*M*, *t*), solid fraction *fαβ*(*M*, *t*) and liquid fraction *fsl*(*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.

As pointed out by [9], in order to ensure the conservativity of the heat flux at the interface

<sup>=</sup> *hi* <sup>+</sup> *hi*<sup>+</sup><sup>1</sup> *hi*<sup>+</sup><sup>1</sup> *λn*+<sup>1</sup> *i*+1

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.

We integrate the differential equation over a space-time finite volume Ω*<sup>i</sup>* × [*tn*, *tn*+1]. A

The backward Euler scheme, first order accurate fully implicit scheme is used. Denoting by

It was proved in [1] that the numerical scheme is *L*<sup>∞</sup> stable, moreover for each time index *n*

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 Ω*i*,*<sup>j</sup>* × [*tn*, *tn*+1]. Each rectangular structured non uniform control volume Ω*i*,*<sup>j</sup>* = [*ri*,*ri*<sup>+</sup>1] × [*zj*, *zj*<sup>+</sup>1],

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

*dz* <sup>=</sup> (*ri*<sup>+</sup>1)

*<sup>g</sup>*+<sup>1</sup> <sup>−</sup> (*ri*)

*<sup>g</sup>* <sup>+</sup> <sup>1</sup> .

*g*+1

*zj*+<sup>1</sup> − *zj*

. (23)

*<sup>i</sup>* the mean value at time *tn* over Ω*<sup>i</sup>* of *ξ*(*x*, *t*), we obtain the following discrete equation

 Ω*i*

*<sup>i</sup>* <sup>+</sup> <sup>Δ</sup>*t kd*(*Tn*+<sup>1</sup>

<sup>1</sup> + <sup>Δ</sup>*t kd*(*Tn*+<sup>1</sup>

*<sup>i</sup>* )

*<sup>i</sup>* )

*<sup>i</sup>* )*<sup>n</sup>* is increasing, hence the discrete finite-volume approximation mimics

 *tn*<sup>+</sup><sup>1</sup> *tn*

+ *hi λn*+<sup>1</sup> *i*

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 507

is computed thanks to the following harmonic

. A decoupled iterative solution of the nonlinear

, (20)

*kd* (*T*) (1 − *ξ*). (21)

*<sup>i</sup>* ∈ [0, 1]. For a given control volume index *i*,

. (22)

*i*+ <sup>1</sup> 2

*λn*+<sup>1</sup> *i*+ <sup>1</sup> 2

*<sup>i</sup>* , *fsl <sup>n</sup>*+<sup>1</sup> *i* 

between two adjacents control volumes, *λn*+<sup>1</sup>

 *Tn*+<sup>1</sup> *<sup>i</sup>* , *<sup>ξ</sup>n*+<sup>1</sup>

cell-centered approximation is used.

*3.1.2. Numerical discretization of the exothermic kinetics*

and for all control volume index *i* = 1, . . . , *I* : *ξ<sup>n</sup>*

**3.2. Two-dimensional finite-volume discretization**

 *zj*<sup>+</sup><sup>1</sup> *zj*

 *ri*<sup>+</sup><sup>1</sup> *ri*

*<sup>r</sup><sup>g</sup> dr*

the behavior of the continuous solution [2].

*d*Ω =

 *tn*<sup>+</sup><sup>1</sup> *tn*

 Ω*i dξ dt* <sup>=</sup>

> *ξn*+<sup>1</sup> *<sup>i</sup>* <sup>=</sup> *<sup>ξ</sup><sup>n</sup>*

mean formula:

where *λn*+<sup>1</sup>

*ξn*

the time sequence (*ξ<sup>n</sup>*

has surface *mi*,*<sup>j</sup>* defined by :

*mi*,*<sup>j</sup>* = Ω*i*,*<sup>j</sup>*

(*g* = 1).

*<sup>i</sup>* = *λ*

#### **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 Ω*<sup>i</sup>* × [*tn*, *tn*+1], where Ω*<sup>i</sup>* = [*ri*,*ri*<sup>+</sup>1] is a control volume. A cell-centered approximation is used. The unknown *T<sup>n</sup> <sup>i</sup>* denotes the mean value of *<sup>T</sup>*(*x*, *<sup>t</sup>*) at time *tn* over <sup>Ω</sup>*i*. *mi* <sup>=</sup> *ri*<sup>+</sup><sup>1</sup> *ri <sup>r</sup><sup>g</sup> dr* is the discrete "mass" of control volume Ω*i*.

$$\int\_{t\_n}^{t\_{n+1}} \int\_{\Omega\_i} \rho \frac{\partial \mathbf{h}\_s \left( T, \tilde{\boldsymbol{\xi}} \right)}{\partial t} \, d\Omega \, dt = \int\_{t\_n}^{t\_{n+1}} \int\_{\Omega\_i} \left( \nabla \left( \lambda \left( T, \tilde{\boldsymbol{\xi}}\_{\boldsymbol{\nu}} f\_{\text{sl}} \right) \nabla \, T \right) + \frac{\rho \Delta H\_f}{M\_{\text{TIC}}} k \left( T \right) \left( 1 - \tilde{\boldsymbol{\xi}} \right) \right) \, d\Omega \, dt. \tag{18}$$

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. Δ*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

$$\begin{split} m\_{i\ell} \rho \frac{h\_{\mathbf{s}}\left(T\_{i}^{n+1}, \mathfrak{z}\_{i}^{n+1}\right) - h\_{\mathbf{s}}\left(T\_{i}^{n}, \mathfrak{z}\_{i}^{n}\right)}{\Delta t} &= \lambda\_{i + \frac{1}{2}}^{n+1} r\_{i+1}^{\mathcal{S}} \frac{T\_{i+1}^{n+1} - T\_{i}^{n+1}}{d\_{i+\frac{1}{2}}} - \lambda\_{i-\frac{1}{2}}^{n+1} r\_{i}^{\mathcal{S}} \frac{T\_{i}^{n+1} - T\_{i-1}^{n+1}}{d\_{i-\frac{1}{2}}} \\ &+ m\_{i} \frac{\rho \Delta H\_{f}}{M\_{\text{TIC}}} k\_{d}(T^{n+1}) \left(1 - \mathfrak{z}\_{i}^{n+1}\right). \end{split} \tag{19}$$

Here *di*<sup>+</sup> <sup>1</sup> 2 is the distance between the center of two adjacent control volumes, and denoting by *hi* = *ri*+<sup>1</sup> − *ri* the length of control volume Ω*i*, then *di*<sup>+</sup> <sup>1</sup> 2 = (*hi*+<sup>1</sup> + *hi*)/2.

As pointed out by [9], in order to ensure the conservativity of the heat flux at the interface between two adjacents control volumes, *λn*+<sup>1</sup> *i*+ <sup>1</sup> 2 is computed thanks to the following harmonic mean formula:

$$
\lambda\_{i+\frac{1}{2}}^{n+1} = \frac{h\_i + h\_{i+1}}{\frac{h\_{i+1}}{\lambda\_{i+1}^{n+1}} + \frac{h\_i}{\lambda\_i^{n+1}}} \,\tag{20}
$$

where *λn*+<sup>1</sup> *<sup>i</sup>* = *λ Tn*+<sup>1</sup> *<sup>i</sup>* , *<sup>ξ</sup>n*+<sup>1</sup> *<sup>i</sup>* , *fsl <sup>n</sup>*+<sup>1</sup> *i* . 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*

6 Numerical Simulations

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

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

We integrate the enthalpy balance over a space-time finite volume Ω*<sup>i</sup>* × [*tn*, *tn*+1], where Ω*<sup>i</sup>* =

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. Δ*t*

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)

> *i*+ <sup>1</sup> 2 *r g i*+1

> > *ρ*Δ*Hf MTiC*

+ *mi*

<sup>∇</sup> (*<sup>λ</sup>* (*T*, *<sup>ξ</sup>*, *fsl*) <sup>∇</sup> *<sup>T</sup>*) <sup>+</sup> *<sup>ρ</sup>*.Δ*Hf*

*Tn*+<sup>1</sup> *<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i di*<sup>+</sup> <sup>1</sup> 2

*kd*(*Tn*+1)

2

is the distance between the center of two adjacent control volumes, and denoting

*<sup>i</sup>* denotes

*d*Ω *dt*. (18)

*ri <sup>r</sup><sup>g</sup> dr* is the discrete "mass" of control

*k* (*T*) (1 − *ξ*)

*Tn*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*−1

*di*<sup>−</sup> <sup>1</sup> 2

. (19)

*MTiC*

<sup>−</sup> *<sup>λ</sup>n*+<sup>1</sup> *<sup>i</sup>*<sup>−</sup> <sup>1</sup> 2 *r g i*

<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>n*+<sup>1</sup> *i* 

= (*hi*+<sup>1</sup> + *hi*)/2.

[*ri*,*ri*<sup>+</sup>1] is a control volume. A cell-centered approximation is used. The unknown *T<sup>n</sup>*

**3. Numerical discretization scheme**

specifics of the enthalpy method.

volume Ω*i*.

 Ω*i ρ*.

is the time-step.

coordinate system

*miρ hs Tn*+<sup>1</sup> *<sup>i</sup>* , *<sup>ξ</sup>n*+<sup>1</sup> *i* − *hs Tn <sup>i</sup>* , *<sup>ξ</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>λ</sup>n*+<sup>1</sup>

Here *di*<sup>+</sup> <sup>1</sup> 2

 *tn*<sup>+</sup><sup>1</sup> *tn*

chemical reaction (differential equation).

*∂hs* (*T*, *ξ*) *∂t*

*3.1.1. Numerical discretization of the enthalpy balance*

the mean value of *<sup>T</sup>*(*x*, *<sup>t</sup>*) at time *tn* over <sup>Ω</sup>*i*. *mi* <sup>=</sup> *ri*<sup>+</sup><sup>1</sup>

 *tn*<sup>+</sup><sup>1</sup> *tn*

 Ω*i* 

*d*Ω *dt* =

by *hi* = *ri*+<sup>1</sup> − *ri* the length of control volume Ω*i*, then *di*<sup>+</sup> <sup>1</sup>

**3.1. One-dimensional finite-volume discretization**

We integrate the differential equation over a space-time finite volume Ω*<sup>i</sup>* × [*tn*, *tn*+1]. A cell-centered approximation is used.

$$\int\_{t\_n}^{t\_{n+1}} \int\_{\Omega\_l} \frac{d\tilde{\xi}}{dt} = \int\_{t\_n}^{t\_{n+1}} \int\_{\Omega\_l} k\_d \left( T \right) \left( 1 - \tilde{\xi} \right) \,. \tag{21}$$

The backward Euler scheme, first order accurate fully implicit scheme is used. Denoting by *ξn <sup>i</sup>* the mean value at time *tn* over Ω*<sup>i</sup>* of *ξ*(*x*, *t*), we obtain the following discrete equation

$$
\tilde{\varsigma}\_{i}^{n+1} = \frac{\tilde{\varsigma}\_{i}^{n} + \Delta t \, k\_{d}(T\_{i}^{n+1})}{1 + \Delta t \, k\_{d}(T\_{i}^{n+1})}.\tag{22}
$$

It was proved in [1] that the numerical scheme is *L*<sup>∞</sup> stable, moreover for each time index *n* and for all control volume index *i* = 1, . . . , *I* : *ξ<sup>n</sup> <sup>i</sup>* ∈ [0, 1]. For a given control volume index *i*, the time sequence (*ξ<sup>n</sup> <sup>i</sup>* )*<sup>n</sup>* 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 Ω*i*,*<sup>j</sup>* × [*tn*, *tn*+1]. Each rectangular structured non uniform control volume Ω*i*,*<sup>j</sup>* = [*ri*,*ri*<sup>+</sup>1] × [*zj*, *zj*<sup>+</sup>1], has surface *mi*,*<sup>j</sup>* defined by :

$$m\_{l,j} = \int\_{\Omega\_{lj}} d\Omega = \int\_{z\_j}^{z\_{j+1}} \left( \int\_{r\_l}^{r\_{l+1}} r^g \, dr \right) \, dz = \frac{(r\_{l+1})^{g+1} - (r\_l)^{g+1}}{g+1} \cdot \left( z\_{j+1} - z\_j \right) . \tag{23}$$

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

8 Numerical Simulations 508 Numerical Simulation – From Theory to Industry Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis <sup>9</sup>

$$\begin{split} & \frac{m\_{i,j}}{\Delta t} \rho. \left( h\_{\sf s} \left( T^{n+1}\_{i,j}, \widetilde{\boldsymbol{\xi}}^{n+1}\_{i,j} \right) - h\_{\sf s} \left( T^{n}\_{\widetilde{i},j}, \widetilde{\boldsymbol{\xi}}^{n}\_{i,j} \right) \right) = \lambda^{n+1}\_{i+\frac{1}{2},j} r^{\mathcal{S}}\_{i+1} \frac{T^{n+1}\_{i+1,j} - T^{n+1}\_{i,j}}{d^{\operatorname{\boldsymbol{\nu}}}\_{i+\frac{1}{2}}} - \lambda^{n+1}\_{i-\frac{1}{2},j} r^{\mathcal{S}}\_{i} \frac{T^{n+1}\_{i,j} - T^{n+1}\_{i-1,j}}{d^{\operatorname{\boldsymbol{\nu}}}\_{i-\frac{1}{2}}} \\ & + \lambda^{n+1}\_{i,j+\frac{1}{2}} \frac{T^{n+1}\_{i,j} - T^{n+1}\_{i,j}}{d^{\operatorname{\boldsymbol{\nu}}}\_{j+\frac{1}{2}}} - \lambda^{n+1}\_{i,j-\frac{1}{2}} \frac{T^{n+1}\_{i,j} - T^{n+1}\_{i,j-1}}{d^{\operatorname{\boldsymbol{\nu}}}\_{j-\frac{1}{2}}} + m\_{i,j} \frac{\rho.\Delta H\_f}{M\_{\operatorname{TdC}}} k \left( T^{n+1}\_{i,j} \right) \left( 1 - \xi^{n+1}\_{i,j} \right). \end{split}$$

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

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 509

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 *Ze* defined

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*∗) =

This subsection presents unpublished results related to the induction time *τind* in

The boundary conditions used to analyze the combustion front propagation in the reactive

We define the induction time, starting point and ending time, three notions that will be heavily

• The induction time *τ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

*T* (0, *t*)

*T* (*Re*, *t*)

4 − *Tfx*<sup>=</sup><sup>0</sup> 4 

4 − *Tfx*<sup>=</sup>*Re* 4 .

*<sup>d</sup>*[*s*−1], *<sup>E</sup>*∗[*kJ*/*mol*]

. *Tad*[*K*] is the adiabatic temperature of the reaction, *Ta*[*K*] is

characterizes the first order exothermic

 <sup>2</sup> <sup>+</sup> <sup>√</sup><sup>5</sup> 

, (25)

, therefore

production.

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

*R T*<sup>2</sup> *ad*

the solid flame propagation is stable.

**4.1. 1D numerical study**

mixture of length *Re* are

used in this section.

*k*∗

(3800 *s*−1, 30.0 *kcal*.*mol*−1), the Zeldovich number *Ze* = 4.52 *< Zec* = 2

planar/cylindrical/spherical geometry such as one depicted in Fig.(21).

*<sup>∂</sup><sup>n</sup>* (0, *<sup>t</sup>*) <sup>=</sup> *� σ*

*<sup>∂</sup><sup>n</sup>* (*Re*, *<sup>t</sup>*) <sup>=</sup> *� σ*

<sup>−</sup>*<sup>λ</sup> <sup>∂</sup><sup>T</sup>*

<sup>−</sup>*<sup>λ</sup> <sup>∂</sup><sup>T</sup>*

*4.1.1. Tools used to analyse the numerical simulation results*

the ambiant temperature and

by *Ze* <sup>=</sup> (*Tad* <sup>−</sup> *Ta*) *<sup>E</sup>*<sup>∗</sup>

The expression used for *λn*+<sup>1</sup> *i*+ <sup>1</sup> 2 ,*j* , *λn*+<sup>1</sup> *<sup>i</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> ,*<sup>j</sup>* , *<sup>λ</sup>n*+<sup>1</sup> *i*,*j*+<sup>1</sup> 2 ,*j* , *λn*+<sup>1</sup> *<sup>i</sup>*,*j*−<sup>1</sup> 2 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 Ω*i*,*j*,*<sup>k</sup>* × [*tn*, *tn*+1]. Each parallelepipedic structured non uniform control volume <sup>Ω</sup>*i*,*j*,*<sup>k</sup>* = [*xi*, *xi*<sup>+</sup>1] <sup>×</sup> [*xj*, *yj*<sup>+</sup>1] <sup>×</sup> [*zk*, *zk*<sup>+</sup>1], has volume *mi*,*j*,*<sup>k</sup>* <sup>=</sup>

Ω*i*,*j*,*<sup>k</sup> d*Ω. *mi*,*j*,*<sup>k</sup>* <sup>Δ</sup>*<sup>t</sup>* .*ρ*. *hs Tn*+<sup>1</sup> *<sup>i</sup>*,*j*,*<sup>k</sup>* , *<sup>ξ</sup>n*+<sup>1</sup> *i*,*j*,*k* − *hs Tn <sup>i</sup>*,*j*,*k*, *<sup>ξ</sup><sup>n</sup> i*,*j*,*k* (24) = *λn*+<sup>1</sup> *i*+ <sup>1</sup> <sup>2</sup> ,*j*,*k Tn*+<sup>1</sup> *<sup>i</sup>*+1,*j*,*<sup>k</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*,*j*,*k dx i*+ <sup>1</sup> 2 <sup>−</sup> *<sup>λ</sup>n*+<sup>1</sup> *<sup>i</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> ,*j*,*k Tn*+<sup>1</sup> *<sup>i</sup>*,*j*,*<sup>k</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*−1,*j*,*k dx <sup>i</sup>*<sup>−</sup> <sup>1</sup> 2 + *λn*+<sup>1</sup> *i*,*j*+<sup>1</sup> 2 ,*k Tn*+<sup>1</sup> *<sup>i</sup>*,*j*+1,*<sup>k</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*,*j*,*k d y j*+ <sup>1</sup> 2 <sup>−</sup> *<sup>λ</sup>n*+<sup>1</sup> *<sup>i</sup>*,*j*<sup>−</sup> <sup>1</sup> 2 ,*k Tn*+<sup>1</sup> *<sup>i</sup>*,*j*,*<sup>k</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*,*j*−1,*k d y <sup>j</sup>*<sup>−</sup> <sup>1</sup> 2 + *λn*+<sup>1</sup> *i*,*j*,*k*+<sup>1</sup> 2 *Tn*+<sup>1</sup> *<sup>i</sup>*,*j*,*k*+<sup>1</sup> <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*,*j*,*k dz k*+ <sup>1</sup> 2 <sup>−</sup> *<sup>λ</sup>n*+<sup>1</sup> *<sup>i</sup>*,*j*,*k*−<sup>1</sup> 2 *Tn*+<sup>1</sup> *<sup>i</sup>*,*j*,*<sup>k</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*,*j*,*k*−1 *dz <sup>k</sup>*<sup>−</sup> <sup>1</sup> 2 + *mi*,*j*,*k*. *ρ*.Δ*Hf MTiC k Tn*+<sup>1</sup> *i*,*j*,*k* <sup>1</sup> <sup>−</sup> *<sup>ξ</sup>n*+<sup>1</sup> *i*,*j*,*k* .

The expression used for *λn*+<sup>1</sup> *i*+ <sup>1</sup> <sup>2</sup> ,*j*,*k* , *λn*+<sup>1</sup> *<sup>i</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> ,*j*,*<sup>k</sup>* , *<sup>λ</sup>n*+<sup>1</sup> *i*,*j*+<sup>1</sup> <sup>2</sup> ,*j*,*k* , *λn*+<sup>1</sup> *<sup>i</sup>*,*j*<sup>−</sup> <sup>1</sup> 2 ,*k* , *λn*+<sup>1</sup> *i*,*j*,*k*+<sup>1</sup> 2 , *λn*+<sup>1</sup> *<sup>i</sup>*,*j*,*k*−<sup>1</sup> 2 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**

8 Numerical Simulations

=*λn*+<sup>1</sup> *i*+ <sup>1</sup> 2 ,*j r g i*+1

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)

We now give the discrete equations related to the finite-volume approximation of the enthalpy balance over a space-time finite volume Ω*i*,*j*,*<sup>k</sup>* × [*tn*, *tn*+1]. Each parallelepipedic structured non uniform control volume <sup>Ω</sup>*i*,*j*,*<sup>k</sup>* = [*xi*, *xi*<sup>+</sup>1] <sup>×</sup> [*xj*, *yj*<sup>+</sup>1] <sup>×</sup> [*zk*, *zk*<sup>+</sup>1], has volume *mi*,*j*,*<sup>k</sup>* <sup>=</sup>

*i*−1,*j*,*k*

*i*,*j*−1,*k*

*i*,*j*,*k*−1

, *λn*+<sup>1</sup> *<sup>i</sup>*,*j*<sup>−</sup> <sup>1</sup> 2 ,*k* , *λn*+<sup>1</sup> *i*,*j*,*k*+<sup>1</sup> 2 , *λn*+<sup>1</sup> *<sup>i</sup>*,*j*,*k*−<sup>1</sup> 2

+ *mi*,*j*,*k*.

*ρ*.Δ*Hf MTiC k Tn*+<sup>1</sup> *i*,*j*,*k* 

+ *mi*,*j*.

*Tn*+<sup>1</sup> *<sup>i</sup>*+1,*<sup>j</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*,*j*

*ρ*.Δ*Hf MTiC k Tn*+<sup>1</sup> *i*,*j*

*dr i*+ <sup>1</sup> 2 <sup>−</sup>*λn*+<sup>1</sup> *<sup>i</sup>*<sup>−</sup> <sup>1</sup> 2 ,*j r g i*

> <sup>1</sup> <sup>−</sup> *<sup>ξ</sup>n*+<sup>1</sup> *i*,*j* .

is a straightforward adaptation of

*Tn*+<sup>1</sup> *<sup>i</sup>*,*<sup>j</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*−1,*j*

*dr <sup>i</sup>*<sup>−</sup> <sup>1</sup> 2

(24)

<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>n*+<sup>1</sup> *i*,*j*,*k* .

is a straightforward

*mi*,*<sup>j</sup>* <sup>Δ</sup>*<sup>t</sup>* .*ρ*. *hs Tn*+<sup>1</sup> *<sup>i</sup>*,*<sup>j</sup>* , *<sup>ξ</sup>n*+<sup>1</sup> *i*,*j* −*hs Tn i*,*j* , *ξ<sup>n</sup> i*,*j* 

+*λn*+<sup>1</sup> *i*,*j*+<sup>1</sup> 2

converges.

Ω*i*,*j*,*<sup>k</sup> d*Ω.

> *mi*,*j*,*<sup>k</sup>* <sup>Δ</sup>*<sup>t</sup>* .*ρ*. *hs Tn*+<sup>1</sup> *<sup>i</sup>*,*j*,*<sup>k</sup>* , *<sup>ξ</sup>n*+<sup>1</sup> *i*,*j*,*k* − *hs Tn <sup>i</sup>*,*j*,*k*, *<sup>ξ</sup><sup>n</sup> i*,*j*,*k*

*Tn*+<sup>1</sup>

*Tn*+<sup>1</sup>

*Tn*+<sup>1</sup>

The expression used for *λn*+<sup>1</sup>

*<sup>i</sup>*+1,*j*,*<sup>k</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*,*j*,*k*

*<sup>i</sup>*,*j*+1,*<sup>k</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*,*j*,*k*

*<sup>i</sup>*,*j*,*k*+<sup>1</sup> <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*,*j*,*k*

*dx i*+ <sup>1</sup> 2

*d y j*+ <sup>1</sup> 2

*dz k*+ <sup>1</sup> 2

**3.4. Numerical implementation**

= *λn*+<sup>1</sup> *i*+ <sup>1</sup> <sup>2</sup> ,*j*,*k*

+ *λn*+<sup>1</sup> *i*,*j*+<sup>1</sup> 2 ,*k*

+ *λn*+<sup>1</sup> *i*,*j*,*k*+<sup>1</sup> 2

*Tn*+<sup>1</sup> *<sup>i</sup>*,*j*+<sup>1</sup> <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*,*j*

> *dz j*+ <sup>1</sup> 2

The expression used for *λn*+<sup>1</sup>

<sup>−</sup> *<sup>λ</sup>n*+<sup>1</sup> *<sup>i</sup>*,*j*−<sup>1</sup> 2

> *i*+ <sup>1</sup> 2 ,*j* , *λn*+<sup>1</sup> *<sup>i</sup>*<sup>−</sup> <sup>1</sup>

**3.3. Three-dimensional finite-volume discretization**

<sup>−</sup> *<sup>λ</sup>n*+<sup>1</sup> *<sup>i</sup>*<sup>−</sup> <sup>1</sup> <sup>2</sup> ,*j*,*k*

<sup>−</sup> *<sup>λ</sup>n*+<sup>1</sup> *<sup>i</sup>*,*j*<sup>−</sup> <sup>1</sup> 2 ,*k*

<sup>−</sup> *<sup>λ</sup>n*+<sup>1</sup> *<sup>i</sup>*,*j*,*k*−<sup>1</sup> 2

> , *λn*+<sup>1</sup> *<sup>i</sup>*<sup>−</sup> <sup>1</sup>

*i*+ <sup>1</sup> <sup>2</sup> ,*j*,*k*

matrix is inverted at each step of the non-linear solver.

*Tn*+<sup>1</sup> *<sup>i</sup>*,*j*,*<sup>k</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup>

*Tn*+<sup>1</sup> *<sup>i</sup>*,*j*,*<sup>k</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup>

*Tn*+<sup>1</sup> *<sup>i</sup>*,*j*,*<sup>k</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup>

<sup>2</sup> ,*j*,*<sup>k</sup>* , *<sup>λ</sup>n*+<sup>1</sup> *i*,*j*+<sup>1</sup> <sup>2</sup> ,*j*,*k*

*dx <sup>i</sup>*<sup>−</sup> <sup>1</sup> 2

*d y <sup>j</sup>*<sup>−</sup> <sup>1</sup> 2

*dz <sup>k</sup>*<sup>−</sup> <sup>1</sup> 2

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

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

*Tn*+<sup>1</sup> *<sup>i</sup>*,*<sup>j</sup>* <sup>−</sup> *<sup>T</sup>n*+<sup>1</sup> *i*,*j*−1

*dz <sup>j</sup>*<sup>−</sup> <sup>1</sup> 2

<sup>2</sup> ,*<sup>j</sup>* , *<sup>λ</sup>n*+<sup>1</sup> *i*,*j*+<sup>1</sup> 2 ,*j* , *λn*+<sup>1</sup> *<sup>i</sup>*,*j*−<sup>1</sup> 2

> 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 *Ze* defined by *Ze* <sup>=</sup> (*Tad* <sup>−</sup> *Ta*) *<sup>E</sup>*<sup>∗</sup> *R T*<sup>2</sup> *ad* . *Tad*[*K*] is the adiabatic temperature of the reaction, *Ta*[*K*] is the ambiant temperature and *k*∗ *<sup>d</sup>*[*s*−1], *<sup>E</sup>*∗[*kJ*/*mol*] 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*∗) = (3800 *s*−1, 30.0 *kcal*.*mol*−1), the Zeldovich number *Ze* = 4.52 *< Zec* = 2 <sup>2</sup> <sup>+</sup> <sup>√</sup><sup>5</sup> , therefore the solid flame propagation is stable.

#### **4.1. 1D numerical study**

This subsection presents unpublished results related to the induction time *τ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 *Re* are

$$\begin{aligned} -\lambda \frac{\partial T}{\partial n}(0, t) &= \varepsilon \,\sigma \left( T \left( 0, t \right)^4 - \left( T\_{f\_{\varepsilon - 0}} \right)^4 \right), \\ -\lambda \frac{\partial T}{\partial n}(R\_{\varepsilon}, t) &= \varepsilon \,\sigma \left( T \left( R\_{\varepsilon}, t \right)^4 - \left( T\_{f\_{\varepsilon - R\_{\varepsilon}}} \right)^4 \right). \end{aligned} \tag{25}$$

#### *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 *τ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* ∈ Ω such that *ξ*(*x*, *τind*) *>* 0. It depends on various parameters such as ignition temperature, heat capacity, thermal conductivity, pre-heating time, mass density.


$$T\_{syn}(\mathbf{x}) = \int\_0^{+\infty} T(\mathbf{x}, t) \frac{d\tilde{\xi}}{dt}(\mathbf{x}, t) \, d\mathbf{x}.\tag{26}$$

**01234567 t [s]**

**01234567 t [s]**

Temperature Field

x [m]

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

0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 <sup>300</sup>

x(ξ=0.01,t) x(ξ=0.50,t) x(ξ=0.99,t) xmax(t)

**Figure 2.** Temporal evolution of *xmax*(*t*) and *x*(*ξ* = 0.01, *t*) *x*(*ξ* = 0.50, *t*) *x*(*ξ* = 0.99, *t*) for wall

**Figure 3.** Temporal evolution of *xmax*(*t*) and *x*(*ξ* = 0.01, *t*) *x*(*ξ* = 0.50, *t*) *x*(*ξ* = 0.99, *t*) for wall

x(ξ=0.01,t) x(ξ=0.50,t) x(ξ=0.99,t) xmax(t)

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 511

**0**

temperature *Tfr*<sup>=</sup><sup>0</sup> = 1600 K. Phase changes are not taken into account.

**0**

**0.002**

temperature *Tfr*<sup>=</sup><sup>0</sup> = 1600 K. Phase changes are taken into account.

T [K]

account.

800

1300

1800

2300

2800

3300

**0.004**

**x [m]**

**0.006**

**0.008**

**0.01**

**0.002**

**0.004**

**x [m]**

**0.006**

**0.008**

**0.01**

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*0.99 *<sup>M</sup>* (*t*) <sup>−</sup> *<sup>x</sup>*0.01 *<sup>M</sup>* (*t*) 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 *xmax*(*t*) gives an information about the thermal history. Fig.(2) and Fig.(3) represent the temporal evolution of *xmax*(*t*) , *x*(*ξ* = 0.01, *t*) *x*(*ξ* = 0.50, *t*) and *x*(*ξ* = 0.99, *t*) for wall temperature *Tfr*<sup>=</sup><sup>0</sup> = 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 *xmax*(*t*) is equal to zero during preheating time, moreover taking into account phase changes as seen on Fig.(3) modifies the velocity of *xmax*(*t*). When the reaction starts, the front propagates, and *xmax*(*t*) evolves differently than *x*0.01 *<sup>M</sup>* (*t*), *<sup>x</sup>*0.99 *<sup>M</sup>* (*t*). Moreover when the front reaches the extremity of the sample, then *xmax*(*t*) 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.

#### *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.

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

#### *4.1.3. Contribution of furnace temperature*

A sample of length *Re* = 1cm is placed inside a furnace maintained at temperature *Tf* ∈ [800, 2400]K. Fig.(6) shows that both induction time *τind* and ending time *τend* increase exponentially when *Tf* decreases, whenever phase change is taken into account or not.

**Figure 2.** Temporal evolution of *xmax*(*t*) and *x*(*ξ* = 0.01, *t*) *x*(*ξ* = 0.50, *t*) *x*(*ξ* = 0.99, *t*) for wall temperature *Tfr*<sup>=</sup><sup>0</sup> = 1600 K. Phase changes are not taken into account.

10 Numerical Simulations

• Thickness of reaction zone. Assuming that the reaction starts at *τind*, then we can follow the evolution of the combustion front propagation in the material thanks to the

• A each instant *t*, *xmax*(*t*) denotes the location which has the maximum temperature in

• The synthesis temperature characterizes the thermal history of the material from the point

A meaningful numerical simulation is presented in Fig.(2) and in Fig.(3). We observe that the

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 *xmax*(*t*) gives an information about the thermal history. Fig.(2) and Fig.(3) represent the temporal evolution of *xmax*(*t*) , *x*(*ξ* = 0.01, *t*) *x*(*ξ* = 0.50, *t*) and *x*(*ξ* = 0.99, *t*) for wall temperature *Tfr*<sup>=</sup><sup>0</sup> = 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 *xmax*(*t*) is equal to zero during preheating time, moreover taking into account phase changes as seen on Fig.(3) modifies the velocity of *xmax*(*t*). When the reaction

when the front reaches the extremity of the sample, then *xmax*(*t*) increases rapidly and soon after decreases because the synthesis is finished and the radiative heat losses contribute

Fig.(4) shows that the propagation of the stiff combustion front is stable, since each profile in

Fig.(5) shows that the energy released by the exothermic kinetics is nearly constant during the

A sample of length *Re* = 1cm is placed inside a furnace maintained at temperature *Tf* ∈ [800, 2400]K. Fig.(6) shows that both induction time *τind* and ending time *τend* increase exponentially when *Tf* decreases, whenever phase change is taken into account or not.

*T*(*x*, *t*)

*<sup>M</sup>* (*t*) <sup>−</sup> *<sup>x</sup>*0.01

*dξ*

*<sup>M</sup>* (*t*), *<sup>t</sup>*) = 0.50 et *<sup>ξ</sup>*(*x*0.99

 +∞ 0

*<sup>M</sup>* (*t*), *<sup>x</sup>*0.50

*<sup>M</sup>* (*t*), *t*) = 0.99,

*<sup>M</sup>* (*t*) et *<sup>x</sup>*0.99

*dt* (*x*, *<sup>t</sup>*) *dx*. (26)

*<sup>M</sup>* (*t*) is nearly constant. Moreover,

*<sup>M</sup>* (*t*), *<sup>x</sup>*0.99

*<sup>M</sup>* (*t*). Moreover

*<sup>M</sup>* (*t*) such

• Starting of the reaction is *xind >* 0 such that *ξ*(*xind*, *τind*) = 0.5,

computation of spatial and temporal evolution of points *x*0.01

*Tsyn*(*x*) =

starts, the front propagates, and *xmax*(*t*) evolves differently than *x*0.01

significantly to the decreasing of the temperature field inside the material.

the steady state propagation regime are equally distant from each other.

propagation of the combustion wave inside the material.

*<sup>M</sup>* (*t*), *<sup>t</sup>*) = 0.01, *<sup>ξ</sup>*(*x*0.50

It was defined and used in [2].

*4.1.2. Combustion front propagation*

*4.1.3. Contribution of furnace temperature*

"thickness" of the reaction-zone defined by *x*0.99

of view of the kinetics and is defined for *x* ∈ Ω by

• Ending time of the reaction is time *τend >* 0 such that ∀*x* ∈ Ω, *ξ*(*x*, *τend*) = 1,

density.

that *ξ*(*x*0.01

domain Ω,

which there exists *x* ∈ Ω such that *ξ*(*x*, *τind*) *>* 0. It depends on various parameters such as ignition temperature, heat capacity, thermal conductivity, pre-heating time, mass

**Figure 3.** Temporal evolution of *xmax*(*t*) and *x*(*ξ* = 0.01, *t*) *x*(*ξ* = 0.50, *t*) *x*(*ξ* = 0.99, *t*) for wall temperature *Tfr*<sup>=</sup><sup>0</sup> = 1600 K. Phase changes are taken into account.

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

**Figure 5.** Temporal evolution of *<sup>ρ</sup>*Δ*Hf MTiC kd*(*T*)(<sup>1</sup> − *<sup>ξ</sup>i*) at equally spaced instants. Phase change is not taken into account.

**0 4 8 12 16 20 24 t [s]**

**0 10 20 30 t [s]**

• 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

Finally in the steady state propagation regime, each curve *xmax*(*t*) 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

Spatial distribution for different values of *Tfr*<sup>=</sup><sup>0</sup> for temperature (*resp. conversion rate*) at *t* = 1s, is represented on Fig.(10). Using a suitable scaling, it is observed that their shape and stiffness

We analyze the time evolution of flux Φ<sup>0</sup> (*t*), exchanged at *r* = 0 between the exterior and the

(0, *<sup>t</sup>*) <sup>=</sup> *ε σ*

can be observed thanks to the analysis of Φ<sup>0</sup> (*t*) represented by Fig.(12) and correlated to the

*<sup>T</sup>*(0, *<sup>t</sup>*)<sup>4</sup> <sup>−</sup> *<sup>T</sup>*<sup>4</sup>

*f* 

. Three consecutive steps

**Tf = 800 K Tf = 1200 K Tf = 1600 K Tf = 2000 K Tf = 2400 K**

**Figure 7.** Temporal evolution of *xind*(*t*) position , for various values of the furnace temperature

**Tf = 800 K Tf = 1200 K Tf = 1600 K Tf = 2000 K Tf = 2400 K**

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 513

**0**

**0.002**

**Figure 8.** Time evolution of *Tmax*(*t*) for several values of temperature furnace

observed that the velocity of the hot spot is a function of *Tf* .

 *λ <sup>∂</sup><sup>T</sup> ∂n* 

**T [K]**

*Tmax*(*t*) depicted in Fig.(8) can be analyzed similarly.

*4.1.4. Energy released/absorbed by the boundary conditions*

**0.004**

**x [m]**

*Tfr*<sup>=</sup><sup>0</sup> ∈ [800, 1200, 1600, 2000, 2400] *K*.

*Tfr*<sup>=</sup><sup>0</sup> ∈ [800, 1200, 1600, 2000, 2400].

material, defined by Φ<sup>0</sup> (*t*) = −

time-evolution of *T*(0, *t*)

is similar.

**0.006**

**0.008**

**0.01**

**Figure 6.** Evolution of *τind* et *τend* as a function of temperature furnace *Tfr*<sup>=</sup><sup>0</sup> .

Fig.(7) shows that the evolution of *xind*(*t*) is nearly similar whatever the furnace temperature *Tf* is. In fact, each curve can be translated from each other, so the behaviour is nearly independent from *Tf* , 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 *Tmax*(*t*). 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.

Fig.(9) shows the same behaviour for *xmax*(*t*) whatever *Tf* is. The following analysis is proposed.


**Figure 7.** Temporal evolution of *xind*(*t*) position , for various values of the furnace temperature *Tfr*<sup>=</sup><sup>0</sup> ∈ [800, 1200, 1600, 2000, 2400] *K*.

12 Numerical Simulations

Energy released by the Kinetics

Energy [J]

0 50000 1E5 1.5E5 2E5 2.5E5 3E5 3.5E5 4E5

**0**

**Figure 6.** Evolution of *τind* et *τend* as a function of temperature furnace *Tfr*<sup>=</sup><sup>0</sup> .

propagation, it reaches the cold extremity, so an intense heat transfer occurs.

constant velocity propagation of TiC synthesis is observed.

**8**

**16**

**time [s]**

**24**

**Figure 5.** Temporal evolution of *<sup>ρ</sup>*Δ*Hf*

into account.

proposed.

is below the ignition temperature.

ending extremity of the sample.

x [m]

τind without phase change τind with phase change τend without phase change τend with phase change

**800 1200 1600 2000 2400 Tf r=0 [K]**

Fig.(7) shows that the evolution of *xind*(*t*) is nearly similar whatever the furnace temperature *Tf* is. In fact, each curve can be translated from each other, so the behaviour is nearly independent from *Tf* , 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 *Tmax*(*t*). 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

Fig.(9) shows the same behaviour for *xmax*(*t*) whatever *Tf* is. The following analysis is

• 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 *xmax*(*t*) = 0, so the temperature

• A combustion front propagation step *xmax*(*t*) increases linearly with respect to time, so a

• An acceleration of the propagation which comes from the influence of the radaitive boundary condition at *x* = *Re*. During a small time interval, the hot spot reaches the

*MTiC kd*(*T*)(<sup>1</sup> − *<sup>ξ</sup>i*) at equally spaced instants. Phase change is not taken

0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02

**Figure 8.** Time evolution of *Tmax*(*t*) for several values of temperature furnace *Tfr*<sup>=</sup><sup>0</sup> ∈ [800, 1200, 1600, 2000, 2400].

• 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 *Tf* .

Finally in the steady state propagation regime, each curve *xmax*(*t*) 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 *Tmax*(*t*) depicted in Fig.(8) can be analyzed similarly.

Spatial distribution for different values of *Tfr*<sup>=</sup><sup>0</sup> for temperature (*resp. conversion rate*) at *t* = 1s, is represented on Fig.(10). Using a suitable scaling, it is observed that their shape and stiffness is similar.

#### *4.1.4. Energy released/absorbed by the boundary conditions*

We analyze the time evolution of flux Φ<sup>0</sup> (*t*), exchanged at *r* = 0 between the exterior and the material, defined by Φ<sup>0</sup> (*t*) = − *λ <sup>∂</sup><sup>T</sup> ∂n* (0, *<sup>t</sup>*) <sup>=</sup> *ε σ <sup>T</sup>*(0, *<sup>t</sup>*)<sup>4</sup> <sup>−</sup> *<sup>T</sup>*<sup>4</sup> *f* . Three consecutive steps can be observed thanks to the analysis of Φ<sup>0</sup> (*t*) represented by Fig.(12) and correlated to the time-evolution of *T*(0, *t*)

**0 10 20 30 t [s]**

**Figure 12.** Time evolution of the energy released/removed by the radiative boundary conditions when

(i) *T*(0, *t*) increases from *Ta* to *Tf* thanks to the heat supplied by radiative heat transfer. Φ<sup>0</sup> (*t*)

(ii) *T*(0, *t*) is above *Tf* and reaches *Tad* (*Tf < Tad*), thanks to the exothermic reaction of titanium-carbide synthesis. The sign of Φ<sup>0</sup> (*t*) changes and its value increases to reach

(iii)*T*(0, *t*) decreases to *Tf* when the kinetics is ended. Φ<sup>0</sup> (*t*) decreases asymptotically to zero.

worth mentioning that due to the thermal shock observed when the "hot" solid combustion front reaches the "cold" boundary, the amplitude Φ*Re* (*t*) is an order of magnitude higher than Φ<sup>0</sup> (*t*). Taking into account phase change doesn't modify the observation done thanks to

We assume that the chemical kinetics has a cut-off temperature *Ts* equal to the first phase transition *Ti<sup>α</sup>* → *Ti<sup>β</sup>* temperature, *Tαβ* = 1166K. We observe on Fig.(13), that the discrepancy

*Tf* <sup>=</sup><sup>900</sup> *<sup>K</sup>* <sup>=</sup>

*end* − *τ*

*ign* increases significantly when *Tf* decreases whenever the phase change is taken

 *τ Ts*=*Tαβ ign*

*Ts*=*Tαβ end* . 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

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

 *λ <sup>∂</sup><sup>T</sup> ∂n* 

**BCr=0(t) : PC = 1 BCr=Re**

**BCr=Re**

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 515

**(t): PC = 0**

**(t): PC = 1**

(*Re*, *<sup>t</sup>*) <sup>=</sup> *ε σ*

*Tf* =1600 *K*

*<sup>T</sup>*(*Re*, *<sup>t</sup>*)<sup>4</sup> <sup>−</sup> *<sup>T</sup>*<sup>4</sup>

. A similar conclusion is

*f* . It is

**2.5e+06 BCr=0(t) : PC = 0**

**−5e+05**

A similar analysis can be done for Φ*Re* (*t*) = −

is negative and its absolute value decreases down to zero,

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

*τTs*=<sup>0</sup> *ign*

phase change is taken into account or not.

a maximum,

Fig.(12).

*τTs*=<sup>0</sup> *ign* − *τ*

*Ts*=*Tαβ*

chemical reaction starts.

into account or not. In practice,

drawn on Fig.(14) for the difference *τTs*=<sup>0</sup>

cut-off temperature is taken into account.

**5e+05**

**1.5e+06**

**Figure 9.** Time evolution of *xmax*(*t*) for different values of furnace temperature *Tfr*<sup>=</sup><sup>0</sup> .

**Figure 10.** Spatial distribution of temperature profile at *t* = 1s, for each wall temperature *Tf* ∈ {800, 1200, 1600, 2000, 2400} K.

**Figure 11.** Spatial distribution of conversion rate profile at *t* = 1s, for each wall temperature *Tf* ∈ {800, 1200, 1600, 2000, 2400} K.

14 Numerical Simulations

**0 10 20 30 t [s]**

> **Tf = 800 K Tf = 1200 K Tf = 1600 K Tf = 2000 K Tf = 2400 K**

**0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 x [m]**

> **Tf = 800 K Tf = 1200 K Tf = 1600 K Tf = 2000 K Tf = 2400 K**

**0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 x [m]**

**Figure 11.** Spatial distribution of conversion rate profile at *t* = 1s, for each wall temperature

**Tf = 800 K Tf = 1200 K Tf = 1600 K Tf = 2000 K Tf = 2400 K**

**0**

**Figure 9.** Time evolution of *xmax*(*t*) for different values of furnace temperature *Tfr*<sup>=</sup><sup>0</sup> .

**Figure 10.** Spatial distribution of temperature profile at *t* = 1s, for each wall temperature

**0.002**

**0**

**0.2**

**0.4**

ξ

**0.6**

**0.8**

**1**

**T [K]**

*Tf* ∈ {800, 1200, 1600, 2000, 2400} K.

*Tf* ∈ {800, 1200, 1600, 2000, 2400} K.

**0.004**

**0.006**

**x [m]**

**0.008**

**0.01**

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


(iii)*T*(0, *t*) decreases to *Tf* when the kinetics is ended. Φ<sup>0</sup> (*t*) decreases asymptotically to zero.

A similar analysis can be done for Φ*Re* (*t*) = − *λ <sup>∂</sup><sup>T</sup> ∂n* (*Re*, *<sup>t</sup>*) <sup>=</sup> *ε σ <sup>T</sup>*(*Re*, *<sup>t</sup>*)<sup>4</sup> <sup>−</sup> *<sup>T</sup>*<sup>4</sup> *f* . It is worth mentioning that due to the thermal shock observed when the "hot" solid combustion front reaches the "cold" boundary, the amplitude Φ*Re* (*t*) is an order of magnitude higher than Φ<sup>0</sup> (*t*). 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 *Ts* equal to the first phase transition *Ti<sup>α</sup>* → *Ti<sup>β</sup>* temperature, *Tαβ* = 1166K. We observe on Fig.(13), that the discrepancy *τTs*=<sup>0</sup> *ign* − *τ Ts*=*Tαβ ign* increases significantly when *Tf* decreases whenever the phase change is taken into account or not. In practice, *τTs*=<sup>0</sup> *ign Tf* <sup>=</sup><sup>900</sup> *<sup>K</sup>* <sup>=</sup> *τ Ts*=*Tαβ ign Tf* =1600 *K* . A similar conclusion is drawn on Fig.(14) for the difference *τTs*=<sup>0</sup> *end* − *τ Ts*=*Tαβ end* .

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.

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.

0 0.0005 0.001 0.0015 0.002 0.0025 0.003 **x [m]**

= 0 K

**Figure 16.** Spatial conversion rate profile at *t* = 2s, for chemical kinetic without cutoff temperature

We consider the situation when we heat identically both extremities of the cylinder with *Tfr*<sup>=</sup><sup>0</sup> = *Tfr*<sup>=</sup>*Re* = *Tf* = 2400 K. Obviously spatial profiles are similar and symmetrical with

> **0 0.002 0.004 0.006 0.008 0.01 x [m]**

Moreover, when the two fronts meet each other at *x* = *Re*/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

Fig.(21) shows the evolution of log(*τind*)(*Tf*) for slab, cylindrical and spherical geometry. The

*ind > τ*

cyl *ind <sup>&</sup>gt; <sup>τ</sup>*sph

*ind* , for each value of

**t= 1 s t= 2 s t= 3 s t= 4 s t= 5 s**

= 1166 K

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 517

ξ(x,t=2 sec) : Ts

ξ(x,t=2 sec) : Ts

0

respect to *x* = *Re*/2 as can be seen on Fig.(17)-(18).

**T [K]**

**Figure 17.** Spatial temperature profiles at several consecutive instants.

case versus the single front case is also noticed in Fig.(20).

three profiles are similar, nevertheless it is observed that *τ*slab

*4.1.7. Contribution of the geometry to the induction time*

*Tf* . Curvature effects may explain this result.

0.2

0.4

0.6

ξ

*Ts* = 0, and with cutoff temperature *Ts* = *Tαβ*. *4.1.6. Analysis of a double front propagation*

0.8

1

**Figure 13.** Evolution of *τind* as a function of *tfr*<sup>=</sup><sup>0</sup> for a kinetics with/without a cut-off temperature. Phase change taken into account (*PC* = 1), or not (*PC* = 0).

**Figure 14.** Evolution of *τind* as a function of *tfr*<sup>=</sup><sup>0</sup> for a kinetics with/without a cut-off temperature. Phase change taken into account (*PC* = 1), or not (*PC* = 0).

**Figure 15.** Spatial temperature profile at *t* = 2s, for chemical kinetic without cutoff temperature *Ts* = 0, and with cutoff temperature *Ts* = *Tαβ*.

**Figure 16.** Spatial conversion rate profile at *t* = 2s, for chemical kinetic without cutoff temperature *Ts* = 0, and with cutoff temperature *Ts* = *Tαβ*.

#### *4.1.6. Analysis of a double front propagation*

16 Numerical Simulations

τign(PC=1,Ts

τign(PC=1,Ts

τign(PC=0,Ts

τign(PC=0,Ts

τend(PC=1,Ts

τend(PC=1,Ts

τend(PC=0,Ts

τend(PC=0,Ts

=0 K)

=0 K)

=1166 K)

=1166 K)

=0 K)

=0 K)

=1166 K)

=1166 K)

**800 1200 1600 2000 2400 Tf r=0 [K]**

**800 1200 1600 2000 2400 Tf r=0 [K]**

0.002 0.0025 0.003 0.0035 0.004 0.0045 0.005 **x [m]**

**Figure 15.** Spatial temperature profile at *t* = 2s, for chemical kinetic without cutoff temperature *Ts* = 0,

**Figure 14.** Evolution of *τind* as a function of *tfr*<sup>=</sup><sup>0</sup> for a kinetics with/without a cut-off temperature.

Ts = 0 K Ts = 1166 K

**Figure 13.** Evolution of *τind* as a function of *tfr*<sup>=</sup><sup>0</sup> for a kinetics with/without a cut-off temperature.

**0**

Phase change taken into account (*PC* = 1), or not (*PC* = 0).

**3**

Phase change taken into account (*PC* = 1), or not (*PC* = 0).

**T [K]**

and with cutoff temperature *Ts* = *Tαβ*.

**9**

**15**

**t [s]**

**21**

**27**

**5**

**10**

**t [s]**

**15**

**20**

**25**

We consider the situation when we heat identically both extremities of the cylinder with *Tfr*<sup>=</sup><sup>0</sup> = *Tfr*<sup>=</sup>*Re* = *Tf* = 2400 K. Obviously spatial profiles are similar and symmetrical with respect to *x* = *Re*/2 as can be seen on Fig.(17)-(18).

**Figure 17.** Spatial temperature profiles at several consecutive instants.

Moreover, when the two fronts meet each other at *x* = *Re*/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).

#### *4.1.7. Contribution of the geometry to the induction time*

Fig.(21) shows the evolution of log(*τind*)(*Tf*) for slab, cylindrical and spherical geometry. The three profiles are similar, nevertheless it is observed that *τ*slab *ind > τ* cyl *ind <sup>&</sup>gt; <sup>τ</sup>*sph *ind* , for each value of *Tf* . Curvature effects may explain this result.

500 1000 1500 2000 2500 **Tf [K]**

**Figure 21.** Evolution of the log(*τind*)(*Tf*) for slab geometry *g* = 0, cylindrical geometry *g* = 1 and

and *kd*(*T*) = *k*<sup>∗</sup>

**Figure 22.** Temperature dependance of *kd*(*T*) for various values of exponent *d* ∈ [0, 2].

*k*∗ <sup>0</sup> = *k*<sup>∗</sup> *d*.*T<sup>d</sup>*

We analyze the contribution of the exponent *d* defined in Eq.(2) to the induction time *τind* and the ending time *τend*. In order to obtain a better precision for *τind*(*d*), we use fractional exponents *d* ∈ [0, 2]. Moreover to be able to compare the results obtained for various values of *kd*(*T*), we perform a normalization of the pre-exponential factor such that its value remain the same when *T* = *Tad* (the adiabatic temperature). We write therefore the equality between

Fig.(22) represents the temperature dependance of *kd*(*T*). Ten numerical simulations are done for five equally-spaced values of *d* ∈ [0, 2]. For each value of *d* we take into account or not the

> **d=0.0 d=0.5 d=1.5 d=1.5 d=2.0**

**300 900 1500 2100 2700 3300 T [K]**

• Each curve *kd*(*T*) is an increasing function of temperature *T* and has the same concavity,

*<sup>d</sup> <sup>T</sup><sup>d</sup> exp*

<sup>−</sup> *<sup>E</sup>*<sup>∗</sup> *R T*

at *T* = *Tad*

*ad*. (27)

g = 0 g = 1 g = 2

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 519

1

*4.1.8. Contribution of the kinetics to the induction time*

<sup>−</sup> *<sup>E</sup>*<sup>∗</sup> *R T*

• When *T* ≤ *Tad*, then *kd*(*T*)*d>*<sup>0</sup> is below *k*0(*T*),

**k(T) [1/s]**

<sup>0</sup> *exp*

phase change and use *Tfr*<sup>=</sup><sup>0</sup> = 1600K.

10

**log (**τign [s])

spherical geometry *g* = 2.

factors *k*0(*T*) = *k*∗

We point out that

**Figure 18.** Spatial conversion rate profiles at several consecutive instants.

**Figure 19.** Temporal evolution of *T*(*Re*, *t*) for single and double fronts propagation.

**Figure 20.** Temporal evolution of *ξ*(*Re*, *t*) for single and double fronts propagation.

**Figure 21.** Evolution of the log(*τind*)(*Tf*) 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*

18 Numerical Simulations

**0 0.002 0.004 0.006 0.008 0.01 x [m]**

**0 1 2 3 4 5 6 7 8 9 10 t [s]**

**0 1 2 3 4 5 6 7 8 9 10 t [s]**

**1 FRONT 2 FRONTS**

**1 FRONT 2 FRONTS**

**t= 1 s t= 2 s t= 3 s t= 4 s t= 5 s**

**0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1**

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

ξ(Re/2, t)

**Figure 19.** Temporal evolution of *T*(*Re*, *t*) for single and double fronts propagation.

**Figure 20.** Temporal evolution of *ξ*(*Re*, *t*) for single and double fronts propagation.

T(Re/2, t) [K]

**Figure 18.** Spatial conversion rate profiles at several consecutive instants.

ξ

We analyze the contribution of the exponent *d* defined in Eq.(2) to the induction time *τind* and the ending time *τend*. In order to obtain a better precision for *τind*(*d*), we use fractional exponents *d* ∈ [0, 2]. Moreover to be able to compare the results obtained for various values of *kd*(*T*), we perform a normalization of the pre-exponential factor such that its value remain the same when *T* = *Tad* (the adiabatic temperature). We write therefore the equality between factors *k*0(*T*) = *k*∗ <sup>0</sup> *exp* <sup>−</sup> *<sup>E</sup>*<sup>∗</sup> *R T* and *kd*(*T*) = *k*<sup>∗</sup> *<sup>d</sup> <sup>T</sup><sup>d</sup> exp* <sup>−</sup> *<sup>E</sup>*<sup>∗</sup> *R T* at *T* = *Tad*

$$k\_0^\* = k\_d^\*. T\_{ad}^d. \tag{27}$$

Fig.(22) represents the temperature dependance of *kd*(*T*). Ten numerical simulations are done for five equally-spaced values of *d* ∈ [0, 2]. For each value of *d* we take into account or not the phase change and use *Tfr*<sup>=</sup><sup>0</sup> = 1600K.

**Figure 22.** Temperature dependance of *kd*(*T*) for various values of exponent *d* ∈ [0, 2]. We point out that


Fig.(23) shows that the increasing rate of the induction time *τind* is super-linear, while it is linear for the ending time *τend* when degree *d* increases. It appears independent from the eventual phase change contribution. It is worth mentioning that all spatial temperature

> **0 0.002 0.004 0.006 0.008 0.01 0.012 x [m]**

**0 0.002 0.004 0.006 0.008 0.01 x [m]**

**Figure 26.** Spatial distribution of synthesis temperature *Tsyn*(*x*,*t* = 30 *s*) for different values of exponent *d*, when phase change is taken into account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K. Fig.(27) shows that the time evolution of *xmax*(*t*) 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 *τind*. Moreover the spatial distribution of the heat released by the kinetics, as seen on Fig.(28) for *t* = 5s presents a maximum for *d* = 0, and a minimum for *d* = 2. This is correlated to the high values of induction time *τind* when high values of *d* are used. Phase change contributes significantly to the time evolution of the front's position *xξ*=1/2(*t*) 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(*t*) is nearly the same. Time evolution of the front's maximum temperature *Tmax*(*t*) 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.

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

d=0.0,PC=0 d=0.0,PC=1 d=0.5,PC=0 d=0.5,PC=1 d=1.0,PC=0 d=1.0,PC=1 d=1.5,PC=0 d=1.5,PC=1 d=2.0,PC=0 d=2.0,PC=1

**Figure 25.** Spatial distribution of conversion rate field *ξ*(*x*, *t* = 5 *s*) for different values of exponent *d*, when phase change is taken into account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K.

d=0.0,PC=0 d=0.0,PC=1 d=0.5,PC=0 d=0.5,PC=1 d=1.0,PC=0 d=1.0,PC=1 d=1.5,PC=0 d=1.5,PC=1 d=2.0,PC=0 d=2.0,PC=1

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 521

**0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1**

**T [K]**

**4.2. 2D numerical study**

ξ

**Figure 23.** Comparison of induction time *τind* and ending time *τend*, when phase change is taken into account (PC=1) or not (PC=0), for several values of exponent *d* ∈ [0, 2], with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K.

profiles represented on Fig.(24) have a discontinuity of the time-derivative when *T* = *Tsl*. This is explained by the contribution *λf us* (*fsl*) of the titanium melting to the thermal conductivity given by Eq.(4).

**Figure 24.** Spatial distribution of temperature field *T*(*x*,*t* = 5 *s*) for different values of exponent *d*, when phase change is taken into account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K.

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

The spatial distribution of the synthesis temperature *Tsyn*(*x*, *t* = 30 *s*) always increases when *d* increases, whether phase change is taken into account or not as seen on Fig.(26).

**Figure 25.** Spatial distribution of conversion rate field *ξ*(*x*, *t* = 5 *s*) for different values of exponent *d*, when phase change is taken into account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K.

**Figure 26.** Spatial distribution of synthesis temperature *Tsyn*(*x*,*t* = 30 *s*) for different values of exponent *d*, when phase change is taken into account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K.

Fig.(27) shows that the time evolution of *xmax*(*t*) 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 *τind*. Moreover the spatial distribution of the heat released by the kinetics, as seen on Fig.(28) for *t* = 5s presents a maximum for *d* = 0, and a minimum for *d* = 2. This is correlated to the high values of induction time *τind* when high values of *d* are used. Phase change contributes significantly to the time evolution of the front's position *xξ*=1/2(*t*) 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(*t*) is nearly the same. Time evolution of the front's maximum temperature *Tmax*(*t*) 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.

#### **4.2. 2D numerical study**

20 Numerical Simulations

• When *T > Tad*, then *kd*(*T*) *>> k*0(*T*). This imply that the heat released by the exothermic kinetics is significantly increasing when *d* increases. Temperature "overshoot" of *Tmax*(*t*)

• We conclude that the velocity of the combustion flame is higher when *d* = 0 than when

Fig.(23) shows that the increasing rate of the induction time *τind* is super-linear, while it is linear for the ending time *τend* when degree *d* increases. It appears independent from the eventual phase change contribution. It is worth mentioning that all spatial temperature

> τind(PC=0) τind(PC=1) τend(PC=0) τend(PC=1)

**0 0.5 1 1.5 2 d**

**Figure 23.** Comparison of induction time *τind* and ending time *τend*, when phase change is taken into account (PC=1) or not (PC=0), for several values of exponent *d* ∈ [0, 2], with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and

profiles represented on Fig.(24) have a discontinuity of the time-derivative when *T* = *Tsl*. This is explained by the contribution *λf us* (*fsl*) of the titanium melting to the thermal conductivity

> **0 0.002 0.004 0.006 0.008 0.01 0.012 x [m]**

**Figure 24.** Spatial distribution of temperature field *T*(*x*,*t* = 5 *s*) for different values of exponent *d*, when phase change is taken into account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K.

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

The spatial distribution of the synthesis temperature *Tsyn*(*x*, *t* = 30 *s*) always increases when

*d* increases, whether phase change is taken into account or not as seen on Fig.(26).

d=0.0,PC=0 d=0.0,PC=1 d=0.5,PC=0 d=0.5,PC=1 d=1.0,PC=0 d=1.0,PC=1 d=1.5,PC=0 d=1.5,PC=1 d=2.0,PC=0 d=2.0,PC=1

*d >* 0, and the temperature obtained is super-adiabatic as seen on Fig.(30).

can be observed on Fig.(30).

*Tfr*<sup>=</sup>*Re* = 300 K.

given by Eq.(4).

account.

**0**

**T [K]**

**5**

**10**

**t [s]**

**15**

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

**0 4 8 12 16 20 t [s]**

**Figure 30.** Time evolution of the front's maximum temperature *Tmax*(*t*) for different values of exponent *d* when phase change is taken into account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K. combustion front inside the cylindrical sample. The boundary conditions are defined over

> � = 0,

�

�

= *ε*.*σ*. �

�

*Tfz*,*r*=*Re*

The computational grid is uniform along *z*-axis and iso-volume along *r*-axis.

(*t*) �

Moreover *T*(*Re*, 0, *t*) *< T*(0, 0, *t*) because of the radiative heat losses at *r* = *Re*.

= *ε*.*σ*. �

= *ε*.*σ*. �

the heat supplied on the exterior surface of the cylinder (radius *Re* = 2cm, height *He* = 2cm) will induce the propagation of the combustion front inside the cylinder, initially at room temperature *T* = *Ta* = 300 *K*. Each computation will analyze the phenomena during *t* = 10 *s*.

thermal diffusion case means that the kinetics is not active, i.e. *kd* = 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 *Ta* ≤ *T*(*r*, *z*, *t*) ≤ *Tf* , 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.

*T* (*Re*, *z*, *t*)

*T* (*r*, *He*, *t*)

4 − �

*T* (*r*, 0, *t*)

(*t*), *Tfr*,*z*=<sup>0</sup> (*t*), *Tfr*,*z*=*He*

4 − � *Tfz*,*r*=*Re* (*t*) �4 � ,

4 − � *Tfr*,*z*=*He* (*t*) �4 � .

> (*t*) �

= (300 *K*, 1600 *K*, 300 *K*) for both simulations. The

*Tfr*,*z*=<sup>0</sup> (*t*)

�4 � ,

characterizes the way

(28)

**d=0.0,PC=0 d=0.0,PC=1 d=0.5,PC=0 d=0.5,PC=1 d=1.0,PC=0 d=1.0,PC=1 d=1.5,PC=0 d=1.5,PC=1 d=2.0,PC=0 d=2.0,PC=1**

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 523

**Tmax(t)**

*∂*Ω = Γ*r*=<sup>0</sup> ∪ Γ*r*=*Re* ∪ Γ*z*=<sup>0</sup> ∪ Γ*z*=*He* by

�

�

�

�

(*t*), *Tfr*,*z*=<sup>0</sup> (*t*), *Tfr*,*z*=*He*

<sup>−</sup>*<sup>λ</sup> <sup>∂</sup><sup>T</sup>* (0, *<sup>z</sup>*, *<sup>t</sup>*) *∂n<sup>e</sup> r*=0

<sup>−</sup>*<sup>λ</sup> <sup>∂</sup><sup>T</sup>* (*r*, 0, *<sup>t</sup>*) *∂n<sup>e</sup> z*=0

<sup>−</sup>*<sup>λ</sup> <sup>∂</sup><sup>T</sup>* (*Re*, *<sup>z</sup>*, *<sup>t</sup>*) *∂n<sup>e</sup> r*=*Re*

<sup>−</sup>*<sup>λ</sup> <sup>∂</sup><sup>T</sup>* (*r*, *He*, *<sup>t</sup>*) *∂n<sup>e</sup> z*=*He*

∀ (*r*, *z*) ∈ Γ*r*=<sup>0</sup> :

∀ (*r*, *z*) ∈ Γ*r*=*Re* :

∀ (*r*, *z*) ∈ Γ*z*=<sup>0</sup> :

∀ (*r*, *z*) ∈ Γ*z*=*He* :

The choice of the temperature triple �

*4.2.1. Reaction-diffusion vs thermal diffusion*

*Tfz*,*r*=*Re*

⎧

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

We choose �

**Figure 27.** Time evolution of *xmax*(*t*) for different values of exponent *d* when phase change is taken into account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 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 *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K.

**Figure 29.** Time evolution of the front's position *xξ*=1/2(*t*) for different values of exponent *d* when phase change is taken into account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K.

**Figure 30.** Time evolution of the front's maximum temperature *Tmax*(*t*) for different values of exponent *d* when phase change is taken into account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K.

combustion front inside the cylindrical sample. The boundary conditions are defined over *∂*Ω = Γ*r*=<sup>0</sup> ∪ Γ*r*=*Re* ∪ Γ*z*=<sup>0</sup> ∪ Γ*z*=*He* by

$$\begin{cases} \forall (r, z) \in \Gamma\_{r=0} \colon \left( -\lambda \frac{\partial \operatorname{T} \left( 0, z, t \right)}{\partial n\_{r}^{e} = 0} \right) = 0, \\ \forall (r, z) \in \Gamma\_{r=R\_{\epsilon}} \colon \left( -\lambda \frac{\partial \operatorname{T} \left( R\_{\epsilon}, z, t \right)}{\partial n\_{r=R\_{\epsilon}}^{e}} \right) = \varepsilon \sigma. \left( \operatorname{T} \left( R\_{\epsilon r} z, t \right)^{4} - \left( T\_{f\_{z|z=\overline{z}\_{\epsilon}}} (t) \right)^{4} \right), \\ \forall (r, z) \in \Gamma\_{z=0} \colon \left( -\lambda \frac{\partial \operatorname{T} \left( r, 0, t \right)}{\partial n\_{z=0}^{e}} \right) = \varepsilon \sigma. \left( \operatorname{T} \left( r, 0, t \right)^{4} - \left( T\_{f\_{z|z=0}} (t) \right)^{4} \right), \\ \forall (r, z) \in \Gamma\_{z=H\_{\epsilon}} \colon \left( -\lambda \frac{\partial \operatorname{T} \left( r, H\_{\epsilon r} t \right)}{\partial n\_{z=H\_{\epsilon}}^{e}} \right) = \varepsilon \sigma. \left( \operatorname{T} \left( r, H\_{\epsilon r} t \right)^{4} - \left( T\_{f\_{z|z=H\_{\epsilon}}} (t) \right)^{4} \right). \end{cases} \tag{28}$$

The choice of the temperature triple � *Tfz*,*r*=*Re* (*t*), *Tfr*,*z*=<sup>0</sup> (*t*), *Tfr*,*z*=*He* (*t*) � characterizes the way the heat supplied on the exterior surface of the cylinder (radius *Re* = 2cm, height *He* = 2cm) will induce the propagation of the combustion front inside the cylinder, initially at room temperature *T* = *Ta* = 300 *K*. Each computation will analyze the phenomena during *t* = 10 *s*. The computational grid is uniform along *z*-axis and iso-volume along *r*-axis.

#### *4.2.1. Reaction-diffusion vs thermal diffusion*

22 Numerical Simulations

d=0.0,PC=0 d=0.0,PC=1 d=0.5,PC=0 d=0.5,PC=1 d=1.0,PC=0 d=1.0,PC=1 d=1.5,PC=0 d=1.5,PC=1 d=2.0,PC=0 d=2.0,PC=1

**1.5 3.5 5.5 7.5 9.5 11.5 13.5 t [s]**

> d=0.0,PC=0 d=0.0,PC=1 d=0.5,PC=0 d=0.5,PC=1 d=1.0,PC=0 d=1.0,PC=1 d=1.5,PC=0 d=1.5,PC=1 d=2.0,PC=0 d=2.0,PC=1

**Figure 27.** Time evolution of *xmax*(*t*) for different values of exponent *d* when phase change is taken into

**0 0.002 0.004 0.006 0.008 0.01 0.012 x [m]**

0 2 4 6 8 10 12 14 **t [s]**

**Figure 29.** Time evolution of the front's position *xξ*=1/2(*t*) for different values of exponent *d* when phase

change is taken into account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 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 *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and

> **d=0.0, PC= 0 d=0.0, PC=1 d=0.5, PC=0 d=0.5, PC=1 d=1.0, PC=0 d=1.0, PC=1 d=1.5, PC=0 d=1.5, PC=1 d=2.0, PC=0 d=2.0, PC=1**

**0**

account (PC=1) or not (PC=0) with *Tfr*<sup>=</sup><sup>0</sup> = 1600 K and *Tfr*<sup>=</sup>*Re* = 300 K.

**0**

**0**

**0.002**

**0.004**

**x [m]**

**0.006**

**0.008**

**0.01**

**2e+05**

**4e+05**

**E**

*Tfr*<sup>=</sup>*Re* = 300 K.

**6e+05**

**8e+05**

**0.002**

**0.004**

**xmax(t) [m]**

**0.006**

**0.008**

**0.01**

We choose � *Tfz*,*r*=*Re* (*t*), *Tfr*,*z*=<sup>0</sup> (*t*), *Tfr*,*z*=*He* (*t*) � = (300 *K*, 1600 *K*, 300 *K*) for both simulations. The thermal diffusion case means that the kinetics is not active, i.e. *kd* = 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 *Ta* ≤ *T*(*r*, *z*, *t*) ≤ *Tf* , 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.

Moreover *T*(*Re*, 0, *t*) *< T*(0, 0, *t*) because of the radiative heat losses at *r* = *Re*.

**Figure 31.** Time evolution of *T*(0, 0, *t*) et *T*(*Re*, 0, *t*).

*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 *Tfz*,*r*=*Re* (*t*), *Tfr*,*z*=<sup>0</sup> (*t*), *Tfr*,*z*=*He* (*t*) = (2400 *K*, 2400 *K*, 300 *K*). 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* = 3s 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 Ω as seen on Fig.(33). A similar conclusion can be drawn upon the shape of the synthesis temperature *Tsyn*(*r*, *z*, 5) represented on Fig.(34), except along the line where the two fronts are joining themselves.

**Figure 33.** Integral of the energy released by the kinetics during the process at each (*r*, *z*) ∈ Ω.

**Figure 34.** Synthesis temperature distribution *Tsyn*(*r*, *z*, 5) for (*r*, *z*) ∈ Ω.

that they can be superposed.

**4.3. 3D numerical study**

at a temperature greater or equal to Θ 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*Θ(*x*) determines if the energy involved in the reaction-diffusion process is uniformly distributed or not in Ω and what is the average temperature of the process. Fig.(35) represents such time distribution at *t* = 5s for Θ = 1800K. Fig.(36) represents such time distribution at *t* = 5s for Θ = 2400K. 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

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 525

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

**Figure 32.** Temperature distribution *T*(*r*, *z*, 3) for (*r*, *z*) ∈ Ω.

According to [2], in order to analyze the kind of propagation, we assume given a temperature Θ ∈ [300, 3000] and a simulation time *tf >* 0. We define for each point *x* ∈ Ω, the thermal history time *t*Θ(*x*) which corresponds to the total time for which a given point *x* ∈ Ω remains

**Figure 33.** Integral of the energy released by the kinetics during the process at each (*r*, *z*) ∈ Ω.

at a temperature greater or equal to Θ 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*Θ(*x*) determines if the energy involved in the reaction-diffusion process is uniformly distributed or not in Ω and what is the average temperature of the process. Fig.(35) represents such time distribution at *t* = 5s for Θ = 1800K. Fig.(36) represents such time distribution at *t* = 5s for Θ = 2400K. 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 *Tsyn*(*r*, *z*, 5) for (*r*, *z*) ∈ Ω.

## **4.3. 3D numerical study**

24 Numerical Simulations

**T(0,0,t) : RD=0,PC=0 T(0,0,t) : RD=0,PC=1 T(0,0,t) : RD=1,PC=0 T(0,0,t) : RD=1,PC=0 T(R,0,t) : RD=0,PC=0 T(R,0,t) : RD=0,PC=1 T(R,0,t) : RD=1,PC=0 T(R,0,t) : RD=1,PC=1**

**0 2 4 6 810 t [s]**

(*t*), *Tfr*,*z*=<sup>0</sup> (*t*), *Tfr*,*z*=*He*

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* = 3s 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 Ω as seen on Fig.(33). A similar conclusion can be drawn upon the shape of the synthesis temperature *Tsyn*(*r*, *z*, 5) represented on Fig.(34), except along the line where the two fronts are joining themselves.

According to [2], in order to analyze the kind of propagation, we assume given a temperature Θ ∈ [300, 3000] and a simulation time *tf >* 0. We define for each point *x* ∈ Ω, the thermal history time *t*Θ(*x*) which corresponds to the total time for which a given point *x* ∈ Ω remains

(*t*) 

= (2400 *K*, 2400 *K*, 300 *K*). At the

*4.2.2. Single, dual longitudinal, dual longitudinal/single radial front* These three cases were considered and analyzed in detail in [2].

> *Tfz*,*r*=*Re*

**Figure 32.** Temperature distribution *T*(*r*, *z*, 3) for (*r*, *z*) ∈ Ω.

**T [K]**

**Figure 31.** Time evolution of *T*(0, 0, *t*) et *T*(*Re*, 0, *t*).

*4.2.3. Dual radial-longitudinal front*

We consider the case where

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

**Figure 35.** Time distribution of the thermal history time *T*1800*K*(*r*, *z*, 5) for (*r*, *z*) ∈ Ω.

**Figure 37.** Temperature distribution *T*(*x*, *y*, *z*, 3) for (*x*, *y*, *z*) ∈ Ω.

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

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis 527

maintains the structured topology of the mesh. It is currently under implementation.

*CEA, Centre DAM-Île de France, Bruyères le Châtel, 91297 Arpajon cedex, France*

*SMS/RMT, LCG, UMR CNRS 5146, Ecole des Mines de Saint-Etienne, 158 cours Fauriel, 42023*

**5. Conclusions and perspectives**

**Author details**

*Saint-Etienne Cedex, France*

A. Aoufi

G. Damamme

**Figure 36.** Time distribution of the thermal history time *T*2400*K*(*r*, *z*, 5) for (*r*, *z*) ∈ Ω.

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*2400*K*(*x*, *y*, *z*, .), *Tsyn*(*x*, *y*, *z*, .).

**Figure 37.** Temperature distribution *T*(*x*, *y*, *z*, 3) for (*x*, *y*, *z*) ∈ Ω.
