Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis

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


Introduction
The purpose of this book chapter is to analyze from a numerical point of view a reaction diffusion mathematical modelling of titanium carbide combustion synthesis from a mixture of titanium and carbide reactive powders thanks to self-propagating high temperature synthesis process.This modelling results in the coupling between a nonlinear parabolic equation expressing the enthalpy balance of the system with radiative boundary conditions and a nonlinear differential equation describing the exothermic chemical reaction in the system.An another multiphysics coupling was analyzed in [3].This Self-propagating High temperature Synthesis (SHS) process was discovered in 1965 by Merzhanov [7], [8] and uses the energy released by the exothermic reaction itself to ensure its self-propagation inside the material after a localized heat supply has been performed for several seconds on the surface of the solid mixture.The stoichiometric solid mixture is made of several kinds of reactive powders.We analyze in this book chapter, the influence of radiative boundary conditions related to the heat supply over the ignition and eventual propagation of a combustion front inside the material.
Four sections are used to present our numerical simulation work.Section two presents the governing equations of the modelling.Section three outline the main aspects of the numerical scheme.Section four analyzes and discusses the main numerical simulation results of the combustion synthesis process.A conclusion summarizes the results that were obtained.

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 ( (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)a t 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.

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 T s = 1166K, which corresponds to the first allotropic phase transition of titanium α,T i α to titanium β,Ti β .An Arrhenius type equation is used ξ(., 0)=0.
The velocity constant k d (T) follows an Arrhenius law in temperature such that In both cases, (with or without cut-off temperature), k * d is a frequency factor in s −1 , E * is the activation energy of the reaction in J/mol, R = 8.314J/(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.

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 •a t T = T α,β = 1166K, the allotropic phase transition of Ti α to Ti β , characterized by its mass latent heat L α,β , [17], and •a t T = T sl = 1943K, titanium melting is characterized by its mass latent heat L sl , [17], and f sl ∈ [0, 1] is the fraction of liquid phase.

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 and the heat capacity of the product C p TiC (T) weighted by the conversion rate of the reaction ξ to obtain where M Ti ,M C 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, ξ, f sl ) Thermal conductivity λ (T, ξ, f sl ) in J.m −1 .K −1 .s−1 is a key parameter governing the propagation of the combustion front.Following [14] we define the effective thermal conductivity λ (T, ξ, f sl ) by where • λ cond (T, ξ) is the thermal conductivity of the system which is a non-linear weighting between the thermal conductivity of the reactants λ Ti+C and the product λ TiC (T a ) according to [17].It is given by while the expression of function f (T) is given in [17].Expression used for λ TiC (T), λ Ti+C (T) are given in [13].
• λ fus ( f sl ) = f sl (1 − ξ) accounts for the appearance of a liquid phase when temperature T = T sl = 1943K, i.e. the temperature for the fusion of titanium.In this case, the fraction f sl 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.This expression represents the fraction of fused titanium that has not yet reacted.
• λ rad (T) expresses the contribution of the radiative heat transfer between the grains and depends on the diameter d p of the titanium/carbide particules, their emissivity ǫ,t h e i r porosity p.The expression used is taken from [18].It was reported in [15] that this term contributes significantly to the velocity of the combustion wave.
• λ conv (T) expresses the contribution of gaz that are present in the porous system.λ air (T) is taken from [13].

Enthalpy balance
The enthalpy balance expressing the local conservation of internal energy can be written as 503 Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis where ρ in kg.m −3 is the density of the system, h the mass enthalpy and λ the thermal conductivity.Temperature T and conversion rate ξ are two thermodynamical variables of the problem.Moreover two scalar variables f αβ and f sl , representing the fraction of solid and liquid phases are also used for the description of the mass enthalpy of the system.We have therefore a set of four independent variables T, ξ, f αβ , f sl that describe the evolution of the system.The mass enthalpy is defined by The computation of the partial derivatives Thanks to the definition of C p (T, ξ) given by Eq.( 3), we obtain where the integral I(T) is computed with the trapezoidal rule.Finally the enthalpy balance for the system is written by where Radiative boundary conditions are defined over ∂Ω for the system and take the form where ǫ = 0.7 is the emissivity of the material while σ is the Stefan-Boltzman constant.The value of T ∞ ∂Ω depends on the boundary ∂Ω.
The initial condition given on Ω states that the sample is at room temperature.Δ r H for T ∈ [300, 3500] K.

Adiabatic temperature of the reaction T ad
The adiabatic temperature T ad of the reaction, or flame temperature, is the maximum temperature when the system is adiabatic.Titanium-carbide is obtained from titanium and carbide elements at temperature T a .The enthalpy of the reaction corresponds to the heat of formation of TiC at temperature T a ,h e n c e(Δ r H) T a = Δ f H 0 (TiC) T a .C o n s i d e r i n g a thermodynamical path composed of the two steps: (i) Titanium and carbide reactants are transformed into titanium-carbide, -TiC-at T a , (ii) End product TiC is heated up from T a to T ad .
One obtains that We use Janaf tables [11] for the expression of the heat capacity as a function of temperature.
When T a = 298K (room temperature value), we solve numerically the integral equation Eq. (17), and obtain that the value of the adiabatic temperature is T ad = 2900K.According to [17], we express

Mathematical modelling of SHS process
The purpose of the numerical simulation is to determine at each point M ∈ Ω,a n da te a c h time t, temperature field T(M, t), conversion rate ξ(M, t), solid fraction f αβ (M, t) and liquid fraction f sl (M, t) that will permit to determine spatial and temporal location of the reaction ignition.The mathematical modelling given by Eq.(1),Eq.(2) and Eq.( 13),Eq.(14),Eq.(15) is well posed and will be solved numerically by the methods described in the next section.

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.

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).

Numerical discretization of the enthalpy balance
We integrate the enthalpy balance over a space-time finite volume ] is a control volume.A cell-centered approximation is used.The unknown T n i denotes the mean value of T(x, t) at time t n over Ω i .m i = r i+1 r i r g dr is the discrete "mass" of control volume Ω i .
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 Here is the distance between the center of two adjacent control volumes, and denoting As pointed out by [9], in order to ensure the conservativity of the heat flux at the interface between two adjacents control volumes, λ n+1 i+ 1   2   is computed thanks to the following harmonic mean formula: where . 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.

Numerical discretization of the exothermic kinetics
We integrate the differential equation over a space-time finite volume The backward Euler scheme, first order accurate fully implicit scheme is used.Denoting by ξ n i the mean value at time t n over Ω i of ξ(x, t), we obtain the following discrete equation It was proved in [1] that the numerical scheme is L ∞ stable, moreover for each time index n and for all control volume index i = 1,...,I : . For a given control volume index i, the time sequence (ξ n i ) n is increasing, hence the discrete finite-volume approximation mimics the behavior of the continuous solution [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 has surface m i,j defined by : The enthalpy balance is written for 2D cartesian geometry (g = 0) or 2D cylindrical geometry (g = 1). .
The expression used for λ n+1 is a straightforward adaptation of Eq.( 20).We use the same decoupled iterative solution of the nonlinear system as in the one-dimensional case.A first-order linearization of both the stiff Arrhenius term, and the enthalpy term is done.This reinforces the strictly dominant sparse penta-diagonal matrix that is inverted at each step of the non-linear solver.It is known that the more strictly dominant a matrix is, the faster iterative solver such as the SSOR, (successive over relaxation method) converges.

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 The expression used for λ n+1 is a straightforward adaptation of Eq.(20).We use the same numerical procedure as the one used for the two-dimensional case.It is worth mentioning that a sparse strictly dominant hepta-diagonal matrix is inverted at each step of the non-linear solver.

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.

1D/2D/3D numerical study
In this section we analyze the ignition and propagation of the combustion front during titanium carbide-TiC-combustion synthesis.It is worth mentioning that the propagation of the combustion wave is either stable i.e. at constant velocity or oscillatory and can be determined according to [16] [8] to sustain this theoretical analysis.Moreover in the constant velocity case resp.(oscillatorycase 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 < Ze c = 2 2 + √ 5 , therefore the solid flame propagation is stable.

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 R e are

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.
• Starting of the reaction is x ind > 0suchthatξ(x ind , τ ind )=0.5, • Ending time of the reaction is time τ end > 0suchthat∀x ∈ Ω, ξ(x, τ end )=1, • Thickness of reaction zone.Assuming that the reaction starts at τ ind ,t h e nw ec a n follow the evolution of the combustion front propagation in the material thanks to the computation of spatial and temporal evolution of points x 0.01 M (t), x 0.50 M (t) et x 0.99 M (t) such that ξ(x 0.01 M (t), t)=0.01,ξ(x 0.50 M (t), t)=0.50 et ξ(x 0.99 M (t), t)=0.99,• A each instant t, x max (t) denotes the location which has the maximum temperature in domain Ω, • The synthesis temperature characterizes the thermal history of the material from the point of view of the kinetics and is defined for x ∈ Ω by 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 M (t) − x 0.01 M (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 x max (t) gives an information about the thermal history.3) modifies the velocity of x max (t).When the reaction starts, the front propagates, and x max (t) evolves differently than x 0.01 M (t), x 0.99 M (t).M o r e o v e r when the front reaches the extremity of the sample, then x max (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.

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.

Contribution of furnace temperature
As a m p l eo fl e n g t hR e = 1cm is placed inside a furnace maintained at temperature T f ∈ [800, 2400]K.Fig. (6) shows that both induction time τ ind and ending time τ end increase exponentially when T f decreases, whenever phase change is taken into account or not.Fig. (7) shows that the evolution of x ind (t) is nearly similar whatever the furnace temperature T f is.In fact, each curve can be translated from each other, so the behaviour is nearly independent from T f , except from small values for which a slight curvature effect is observed.
A similar conclusion can be drawn upon analyzing Fig. (8) which represents the time evolution of T max (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 x max(t) whatever T f is.The following analysis is proposed.
• A pre-heating stage, for which thanks to the heat supply at x = 0 (radiative boundary condition), there is a gradual increase in temperature.But x max (t)=0, so the temperature is below the ignition temperature.
• A combustion front propagation step x max (t) increases linearly with respect to time, so a constant velocity propagation of TiC synthesis is observed.
• An acceleration of the propagation which comes from the influence of the radaitive boundary condition at x = R e .During a small time interval, the hot spot reaches the ending extremity of the sample.• A cooling stage by thermal conduction, since the combustion synthesis is finished.The hot spot moves rapidly towards the center of the material thanks to the cooling.It is also observed that the velocity of the hot spot is a function of T f .
Finally in the steady state propagation regime, each curve x max (t) can be translated from each other.The local heat supply induced by the boundary condition at x = 0doesn'tmodifythe characteristics of the combustion front propagation such as it's velocity.The time evolution of T max(t) depicted in Fig.( 8) can be analyzed similarly.
Spatial distribution for different values of T f r=0 for temperature (resp.conversion rate)att = 1s, is represented on Fig. (10).Using a suitable scaling, it is observed that their shape and stiffness is similar.

Energy released/absorbed by the boundary conditions
We analyze the time evolution of flux Φ 0 (t), exchanged at r = 0 between the exterior and the material, defined by Φ 0 (t) = − λ ∂T ∂n (0, t) = εσ T(0, t) 4 − T 4 f .Three consecutive steps can be observed thanks to the analysis of Φ 0 (t) represented by Fig. (12) and correlated to the time-evolution of T(0, t)  (i) T(0, t) increases from T a to T f thanks to the heat supplied by radiative heat transfer.Φ 0 (t) is negative and its absolute value decreases down to zero, (ii) T(0, t) is above T f and reaches T ad (T f < T ad ), thanks to the exothermic reaction of titanium-carbide synthesis.The sign of Φ 0 (t) changes and its value increases to reach amaximum, (iii)T(0, t) decreases to T f when the kinetics is ended.Φ 0 (t) decreases asymptotically to zero.
A similar analysis can be done for Φ R e (t) = − λ ∂T ∂n (R e , t) = εσ T(R e , t) 4 − T 4 f .I ti s worth mentioning that due to the thermal shock observed when the "hot" solid combustion front reaches the "cold" boundary, the amplitude Φ R e (t) is an order of magnitude higher than Φ 0 (t).Taking into account phase change doesn't modify the observation done thanks to Fig. (12).

Contribution of the cut-off temperature to the ignition and propagation
We assume that the chemical kinetics has a cut-off temperature T s equal to the first phase transition Ti α → Ti β temperature, T αβ = 1166K.We observe on Fig. (13), that the discrepancy increases significantly when T f decreases whenever the phase change is taken into account or not.In practice, 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.τ ign (PC=1,T s =0 K) τ ign (PC=1,T s =1166 K) τ ign (PC=0,T s =0 K) τ ign (PC=0,T s =1166 K)

Analysis of a double front propagation
We consider the situation when we heat identically both extremities of the cylinder with T f r=0 = T f r=Re = T f = 2400 K. Obviously spatial profiles are similar and symmetrical with respect to x = R e /2 as can be seen on Fig.( 17)-( 18).

Contribution of the kinetics to the induction time
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 k d (T), we perform a normalization of the pre-exponential factor such that its value remain the same when T = T ad (the adiabatic temperature).We write therefore the equality between We point out that •E a c h c u r v e k d (T) is an increasing function of temperature T and has the same concavity, •W h e n T ≤ T ad ,thenk d (T) d>0 is below k 0 (T), •W h e n T > T ad ,thenk d (T) >> k 0 (T).This imply that the heat released by the exothermic kinetics is significantly increasing when d increases.Temperature "overshoot" of T max (t) can be observed on Fig.

(30).
• We conclude that the velocity of the combustion flame is higher when d = 0 than when d > 0, and the temperature obtained is super-adiabatic as seen on Fig. (30).30) is in fact nearly independent from d when a suitable time-translation is performed whenever phase change is taken into account or not.

2D numerical study
This subsection computes the propagation of a dual radial/longitudinal front and a single longitudinal front for various values of the temperature furnace in which the cylindrical sample is placed.The main interest of this modelling with respect to the 1D case is to analyze the contribution of the lateral heat losses over the ignition and propagation of the   combustion front inside the cylindrical sample.The boundary conditions are defined over The choice of the temperature triple T f z,r=Re (t), T f r,z=0 (t), T f r,z= He (t) characterizes the way the heat supplied on the exterior surface of the cylinder (radius R e = 2cm, height H e = 2cm) will induce the propagation of the combustion front inside the cylinder, initially at room temperature T = T a = 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.

Reaction-diffusion vs thermal diffusion
We choose T f z,r=Re (t), T f r,z=0 (t), T f r,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. k d = 0, so the initial mixture of reactive powders is not transformed.At each point (r, z) of the sample and at each instant t > 0, we observe on Fig.( 31) that T a ≤ T(r, z, t) ≤ T f , when the exothermic kinetics is not taken into account, moreover the numerical solution fulfills a discrete maximum principle.A similar principle is observed when the kinetics is taken into account, but the time evolution is significantly different because of the sharp rise in temperature when the synthesis starts.
Moreover T(R e ,0,t) < T(0, 0, t) because of the radiative heat losses at r = R e .

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

Dual radial-longitudinal front
We consider the case where T f z,r=Re (t), T f r,z=0 (t), T f r,z= He (t) = (2400 K, 2400 K, 300 K).A tt h e 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  According to [2], in order to analyze the kind of propagation, we assume given a temperature Θ ∈ [300, 3000] and a simulation time t f > 0. We define for each point x ∈ Ω,t h et h e r m a l history time t Θ (x) which corresponds to the total time for which a given point x ∈ Ω remains

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

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

Fig. 1
Fig.1 represents the evolution of I(T) Δ r H for T ∈ [300, 3500].I t i s o b s e r v e d t h a t

Figure 5 .τFigure 6 .
Figure 5. Temporal evolution of ρΔH f M TiC k d (T) (1 − ξ i ) at equally spaced instants.Phase change is not taken into account.

Figure 9 .Figure 10 .Figure 11 .
Figure 9.Time evolution of x max (t) for different values of furnace temperature T f r=0 .

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

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

Figure 14 .Figure 15 .Figure 16 .
Figure 14.Evolution of τ ind as a function of t f r=0 for a kinetics with/without a cut-off temperature.Phase change taken into account (PC = 1), or not (PC = 0).

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

Fig.( 21
Fig.(21) shows the evolution of log(τ ind )(T f ) for slab, cylindrical and spherical geometry.The three profiles are similar, nevertheless it is observed that τ slab ind > τ cyl ind > τ sph ind ,foreachvalueof T f .Curvature effects may explain this result.

Fig.( 22 Figure 22 .
Fig.(22) represents the temperature dependance of k d (T).Ten numerical simulations are done for five equally-spaced values of d ∈ [0, 2].Forea c hv a lueofd we take into account or not the phase change and use T f r=0 = 1600K.

Fig.( 23 Figure 23 .
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

Figure 24 .
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 T f r=0 = 1600 K and T f r=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 T syn (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 .Figure 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 T f r=0 = 1600 K and T f r=Re = 300 K.

Fig.( 27
Fig.(27) shows that the time evolution of x max (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 T max (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.

521Figure 27 .Figure 28 .
Figure 27.Time evolution of x max (t) for different values of exponent d when phase change is taken into account (PC=1) or not (PC=0) with T f r=0 = 1600 K and T f r=Re = 300 K.

Figure 29 .
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 T f r=0 = 1600 K and T f r=Re = 300 K.

Figure 30 .
Figure 30.Time evolution of the front's maximum temperature T max (t) for different values of exponent d when phase change is taken into account (PC=1) or not (PC=0) with T f r=0 = 1600 K and T f r=Re = 300 K.
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 T syn (r, z,5) represented on Fig.(34), except along the line where the two fronts are joining themselves.

Figure 33 .
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 4.16s (resp.4.70) for 1800K (resp.2400K), and that they can be superposed.
by the computation of the Zeldovich number Ze defined ad − T a ) E * RT 2 ad .T ad [K] is the adiabatic temperature of the reaction, T a [K] is the ambiant temperature and k * d [s −1 ], E * [kJ/mol] characterizes the first order exothermic kinetics.R is the gaz constant.Experimental observations were reported in Temporal evolution of x max (t) and x(ξ = 0.01, t) x(ξ = 0.50, t) x(ξ = 0.99, t) for wall temperature T f r=0 = 1600 K. Phase changes are not taken into account.