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 . This Self-propagating High temperature Synthesis (SHS) process was discovered in 1965 by Merzhanov ,  and uses the energy released by the exothermic reaction itself to ensure its self-propagation inside the material after a localized heat supply has been performed for several seconds on the surface of the solid mixture. The stoichiometric solid mixture is made of several kinds of reactive powders. We analyze in this book chapter, the influence of radiative boundary conditions related to the heat supply over the ignition and eventual propagation of a combustion front inside the material.
Four sections are used to present our numerical simulation work. Section two presents the governing equations of the modelling. Section three outline the main aspects of the numerical scheme. Section four analyzes and discusses the main numerical simulation results of the combustion synthesis process. A conclusion summarizes the results that were obtained.
2. Mathematical modelling
This section describes the main features of the modelling which expresses the coupling between a reaction-diffusion written for the enthalpy balance and a differential equation written for the exothermic chemical kinetics. SHS (Self-propagation High-temperature Synthesis) process is a condensed phase process which converts a mixture of powders into an end product. In this paper we consider the synthesis of titanium carbide in solid phase ( subscript) thanks to the exothermic reaction . We assume that we have a stoichiometric and isotropic mixture of compacted powders of titanium and carbide. Let us denote by ( resp. ) the temperature ( resp. conversion rate) at point at time . The system is composed of the fraction of titanium carbide TiC end product and the fraction of remaining powders of titanium and carbide.
A first order exothermic kinetics is used to describe the synthesis of titanium carbide, hence a single variable, the conversion rate of the reaction is defined. When ( resp. ) the reaction has not started ( resp. is ended). The velocity of the reaction may have or not a cutoff temperature K, which corresponds to the first allotropic phase transition of titanium , Tito titanium , Ti. An Arrhenius type equation is used
The velocity constant follows an Arrhenius law in temperature such that
In both cases, (with or without cut-off temperature), is a frequency factor in , is the activation energy of the reaction in , is the perfect gas constant and 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 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 and the heat capacity of the product weighted by the conversion rate of the reaction to obtain
2.2.2. Equation used for the thermal conductivity
Thermal conductivity in is a key parameter governing the propagation of the combustion front. Following  we define the effective thermal conductivity by
2.2.3. Enthalpy balance
The enthalpy balance expressing the local conservation of internal energy can be written as
where in is the density of the system, the mass enthalpy and the thermal conductivity. Temperature and conversion rate are two thermodynamical variables of the problem. Moreover two scalar variables and , 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 that describe the evolution of the system. The mass enthalpy is defined by
The computation of the partial derivatives , , and leads to
Thanks to the definition of given by Eq.(3), we obtain
where the integral is computed with the trapezoidal rule.
Fig.1 represents the evolution of for . It is observed that . The contribution of to is therefore neglected.
Finally the enthalpy balance for the system is written by
Radiative boundary conditions are defined over for the system and take the form
where is the emissivity of the material while is the Stefan-Boltzman constant. The value of depends on the boundary .
The initial condition given on states that the sample is at room temperature.
2.2.4. Adiabatic temperature of the reaction
The adiabatic temperature 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 . The enthalpy of the reaction corresponds to the heat of formation of TiC at temperature , hence . Considering a thermodynamical path composed of the two steps:
Titanium and carbide reactants are transformed into titanium-carbide, -TiC- at ,
End product TiC is heated up from to .
One obtains that
We use Janaf tables  for the expression of the heat capacity as a function of temperature. When (room temperature value), we solve numerically the integral equation Eq.(17), and obtain that the value of the adiabatic temperature is . According to , we express .
2.3. Mathematical modelling of SHS process
The purpose of the numerical simulation is to determine at each point , and at each time , temperature field , conversion rate , solid fraction and liquid fraction 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 . Error estimates are given in . 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 . This will be omitted from now, since we will mainly focuss our attention on the construction of the numerical scheme for the reaction-diffusion equation on a non uniform structured mesh and discard formally the specifics of the enthalpy method.
3.1. One-dimensional finite-volume discretization
We present the set of discrete equations that arise from the implicit finite-volume discretization of the enthalpy balance (reaction-diffusion equation) and the exothermic chemical reaction (differential equation).
3.1.1. Numerical discretization of the enthalpy balance
We integrate the enthalpy balance over a space-time finite volume , where is a control volume. A cell-centered approximation is used. The unknown denotes the mean value of at time over . is the discrete "mass" of control volume .
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. is the time-step.
We give only for internal control volumes, the finite-volume discretization of the enthalpy balance in a one-dimensional cartesian (), cylindrical () and spherical () coordinate system
Here is the distance between the center of two adjacent control volumes, and denoting by the length of control volume , then .
As pointed out by , in order to ensure the conservativity of the heat flux at the interface between two adjacents control volumes, 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 . The numerical solution cost of the three-diagonal linear system by this method increases linearly with the number of unknowns.
3.1.2. Numerical discretization of the exothermic kinetics
We integrate the differential equation over a space-time finite volume . A cell-centered approximation is used.
The backward Euler scheme, first order accurate fully implicit scheme is used. Denoting by the mean value at time over of , we obtain the following discrete equation
It was proved in  that the numerical scheme is stable, moreover for each time index and for all control volume index . For a given control volume index , the time sequence is increasing, hence the discrete finite-volume approximation mimics the behavior of the continuous solution .
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 . Each rectangular structured non uniform control volume , has surface defined by :
The enthalpy balance is written for 2D cartesian geometry or 2D cylindrical geometry .
The expression used for , , , 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 . Each parallelepipedic structured non uniform control volume , has volume .
The expression used for , , , , , 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 . 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 . Efforts were done to achieve the highest possible performance on processors that use memory cache hierarchy such as MIPS R1X00 processors and INTEL Core i7 processors. As an example, all loops invariant quantities are pre-computed and stored into one-dimensional arrays. Dynamically allocated multidimensional arrays are not used, because of pointer aliasing problems. We prefer to convert such arrays into one-dimensional arrays. A storage similar to CRS(Compact Row Storage) is used that takes into account the penta (hepta)-diagonal sparsity of the matrices that are iteratively inverted by SSOR and SIP methods. The postprocessing of the software’s results saved to vtk format is done thanks to mayavi software for temporal snapshot and Paraview for interactive analysis and movie production.
4. 1D/2D/3D numerical study
In this section we analyze the ignition and propagation of the combustion front during titanium carbide-TiC- combustion synthesis. It is worth mentioning that the propagation of the combustion wave is either stable i.e. at constant velocity or oscillatory and can be determined according to  by the computation of the Zeldovich number defined by . is the adiabatic temperature of the reaction, is the ambiant temperature and characterizes the first order exothermic kinetics. is the gaz constant. Experimental observations were reported in  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 , the Zeldovich number , therefore the solid flame propagation is stable.
4.1. 1D numerical study
This subsection presents unpublished results related to the induction time 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 are
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 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 such that . It depends on various parameters such as ignition temperature, heat capacity, thermal conductivity, pre-heating time, mass density.
Starting of the reaction is such that ,
Ending time of the reaction is time such that ,
Thickness of reaction zone. Assuming that the reaction starts at , then we can follow the evolution of the combustion front propagation in the material thanks to the computation of spatial and temporal evolution of points , et such that , et ,
A each instant , 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 by
It was defined and used in .
A meaningful numerical simulation is presented in Fig.(2) and in Fig.(3). We observe that the "thickness" of the reaction-zone defined by 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 gives an information about the thermal history. Fig.(2) and Fig.(3) represent the temporal evolution of , and for wall temperature 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 is equal to zero during preheating time, moreover taking into account phase changes as seen on Fig.(3) modifies the velocity of . When the reaction starts, the front propagates, and evolves differently than . Moreover when the front reaches the extremity of the sample, then 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 cm is placed inside a furnace maintained at temperature K. Fig.(6) shows that both induction time and ending time increase exponentially when decreases, whenever phase change is taken into account or not.
Fig.(7) shows that the evolution of is nearly similar whatever the furnace temperature is. In fact, each curve can be translated from each other, so the behaviour is nearly independent from , 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 . 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 whatever is. The following analysis is proposed.
A pre-heating stage, for which thanks to the heat supply at (radiative boundary condition), there is a gradual increase in temperature. But , so the temperature is below the ignition temperature.
A combustion front propagation step 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 . 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 .
Finally in the steady state propagation regime, each curve can be translated from each other. The local heat supply induced by the boundary condition at doesn’t modify the characteristics of the combustion front propagation such as it’s velocity. The time evolution of depicted in Fig.(8) can be analyzed similarly.
Spatial distribution for different values of for temperature ( resp. conversion rate) at s, 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 , exchanged at between the exterior and the material, defined by . Three consecutive steps can be observed thanks to the analysis of represented by Fig.(12) and correlated to the time-evolution of
A similar analysis can be done for . It is worth mentioning that due to the thermal shock observed when the "hot" solid combustion front reaches the "cold" boundary, the amplitude is an order of magnitude higher than . 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 equal to the first phase transition temperature, K. We observe on Fig.(13), that the discrepancy increases significantly when decreases whenever the phase change is taken into account or not. In practice, . A similar conclusion is drawn on Fig.(14) for the difference .
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.
4.1.6. Analysis of a double front propagation
We consider the situation when we heat identically both extremities of the cylinder with K. Obviously spatial profiles are similar and symmetrical with respect to as can be seen on Fig.(17)-(18).
Moreover, when the two fronts meet each other at , 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 for slab, cylindrical and spherical geometry. The three profiles are similar, nevertheless it is observed that , for each value of . Curvature effects may explain this result.
4.1.8. Contribution of the kinetics to the induction time
We analyze the contribution of the exponent defined in Eq.(2) to the induction time and the ending time . In order to obtain a better precision for , we use fractional exponents . Moreover to be able to compare the results obtained for various values of , we perform a normalization of the pre-exponential factor such that its value remain the same when (the adiabatic temperature). We write therefore the equality between factors and at
Fig.(22) represents the temperature dependance of . Ten numerical simulations are done for five equally-spaced values of . For each value of we take into account or not the phase change and use K.
We point out that
Each curve is an increasing function of temperature and has the same concavity,
When , then is below ,
When , then . This imply that the heat released by the exothermic kinetics is significantly increasing when increases. Temperature "overshoot" of can be observed on Fig.(30).
We conclude that the velocity of the combustion flame is higher when than when , and the temperature obtained is super-adiabatic as seen on Fig.(30).
Fig.(23) shows that the increasing rate of the induction time is super-linear, while it is linear for the ending time when degree increases. It appears independent from the eventual phase change contribution.
It is worth mentioning that all spatial temperature profiles represented on Fig.(24) have a discontinuity of the time-derivative when . This is explained by the contribution of the titanium melting to the thermal conductivity given by Eq.(4).
Conversion rate profiles represented on Fig.(25) are sharper when phase change is taken into account.
The spatial distribution of the synthesis temperature always increases when increases, whether phase change is taken into account or not as seen on Fig.(26).
Fig.(27) shows that the time evolution of has a similar shape for each value of phase change is taken into account or not. Considering high values of imply a significant delay for . Moreover the spatial distribution of the heat released by the kinetics, as seen on Fig.(28) for s presents a maximum for , and a minimum for . This is correlated to the high values of induction time when high values of are used. Phase change contributes significantly to the time evolution of the front’s position as observed on Fig.(29). Without phase change, the evolution is nearly independent from the value of since the slope of is nearly the same. Time evolution of the front’s maximum temperature for different values of exponent as seen on Fig.(30) is in fact nearly independent from when a suitable time-translation is performed whenever phase change is taken into account or not.
4.2. 2D numerical study
This subsection computes the propagation of a dual radial/longitudinal front and a single longitudinal front for various values of the temperature furnace in which the cylindrical sample is placed. The main interest of this modelling with respect to the 1D case is to analyze the contribution of the lateral heat losses over the ignition and propagation of the combustion front inside the cylindrical sample. The boundary conditions are defined over by
The choice of the temperature triple characterizes the way the heat supplied on the exterior surface of the cylinder (radius cm, height cm) will induce the propagation of the combustion front inside the cylinder, initially at room temperature . Each computation will analyze the phenomena during . The computational grid is uniform along -axis and iso-volume along -axis.
4.2.1. Reaction-diffusion vs thermal diffusion
We choose for both simulations. The thermal diffusion case means that the kinetics is not active, i.e. , so the initial mixture of reactive powders is not transformed. At each point of the sample and at each instant , we observe on Fig.(31) that , 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 because of the radiative heat losses at .
4.2.2. Single, dual longitudinal, dual longitudinal/single radial front
These three cases were considered and analyzed in detail in .
4.2.3. Dual radial-longitudinal front
We consider the case where . 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 s represented by Fig.(32). The conversion rate 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 represented on Fig.(34), except along the line where the two fronts are joining themselves.
According to , in order to analyze the kind of propagation, we assume given a temperature and a simulation time . We define for each point , the thermal history time which corresponds to the total time for which a given point remains 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, 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 s for K. Fig.(36) represents such time distribution at s for K. As in the previous figures, it is influenced by the propagation of the combustion wave. It is worth mentioning that on these two figures the time distribution is nearly equal, from up to 4.16s (resp. 4.70) for 1800K (resp. 2400K), and that they can be superposed.
4.3. 3D numerical study
This subsection accounts for the heterogeneity of the radiative boundary conditions to analyze the ignition and propagation of the combustion front in a cubic sample. We assume that the heat supplied to the six faces of the cubic sample is different, in the sense that the wall temperature that contributes to the preheating of each face of the cubic sample is different. This will induce a specific transient spatial pattern of the combustion front. A meaningful snapshot of the transient temperature profile is depicted in Fig.(37). It is clearly observed that the shape of the front’s propagation inside the cube is influenced by the boundary conditions. If a regular propagation is required, providing heat supply over one single face of the cube is enough to obtain quickly titanium-carbide. So combining a different heat supply system for each face of the cube appears complicated from an experimental point of view. Numerical simulation show that the pattern of the propagation in this case is more complicated than in the single heat supply case and requires more computational power to get finely resolved. An analysis, not presented here due space constraints, and using the same methodology as in the previous 2D case show that the same conclusions can be drawn for
5. Conclusions and perspectives
In this book chapter, we have presented a multidimensional modelling for the numerical computation of ignition and propagation of combustion fronts during combustion synthesis of ceramic materials such as titanium-carbide. A detailed computational study was done in 1D slab/cylindrical/spherical geometry, 2D cylindrical/cartesian and 3D cartesian geometry to analyze the influence of the radiative boundary conditions over the induction time. The radiation contribution to the thermal conductivity was taken into account and the sensitivity of the induction time to several parameters such as the kinetics, wall temperature, phase change was carefully analyzed. Our numerical software Hephaïstos was presented. It implements an implicit finite-volume scheme for which error estimates, discrete maximum principles were reported and used to ensure the consistency of the numerical results. This modelling study was done for titanium carbide TiC. It can be applied to other ceramic materials such as silicium carbide-SiC-. Moreover in order to analyze the combustion front propagation for high Zeldovich number cases, an adaptive finite-volume scheme is required. A new monotonicity preserving refinement/derefinement conservative algorithm has been designed for multidimensional computations in various coordinate systems. This algorithm maintains the structured topology of the mesh. It is currently under implementation.