Open access peer-reviewed chapter

Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis

By A. Aoufi and G. Damamme

Submitted: November 28th 2011Reviewed: May 3rd 2012Published: September 19th 2012

DOI: 10.5772/48441

Downloaded: 1085

1. Introduction

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

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

2. Mathematical modelling

This section describes the main features of the modelling which expresses the coupling between a reaction-diffusion written for the enthalpy balance and a differential equation written for the exothermic chemical kinetics. SHS (Self-propagation High-temperature Synthesis) process is a condensed phase process which converts a mixture of powders into an end product. In this paper we consider the synthesis of titanium carbide in solid phase ( (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=TM,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.

A first order exothermic kinetics is used to describe the synthesis of titanium carbide, hence a single variable, the conversion rate Ϝ0,1of 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Ϗto titanium ϐ, Tiϐ. An Arrhenius type equation is used

dϜdt=kdT1-Ϝ,Ϝ(.,0)=0.E1

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

kdT=kd*Tdexp-E*R.T.E2

In both cases, (with or without cut-off temperature), kd*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 d0,2is 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=4000Kis high and gives an idea of the stiffness of the reaction-diffusion/kinetics coupling.

2.2. Heat transfer modelling

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

  • T=TϏ,ϐ=1166TiϏTiϐLϏ,ϐ17fϏϐ[0,1]ϐ
  • T=Tsl=1943Lsl17fsl[0,1]

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 CpTiT,CpCTand the heat capacity of the product CpTiCTweighted by the conversion rate of the reaction Ϝto obtain

CpT,Ϝ=Ϝ.CpTiCT+1-Ϝ.MTi.CpTiT+MC.CpCTMTiC,E3

where MTi,MCet MTiCdenote 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,Ϝ,fslin J.m-1.K-1.s-1is a key parameter governing the propagation of the combustion front. Following [14] we define the effective thermal conductivity ϙT,Ϝ,fslby

ϙT,Ϝ,fsl=ϙcondT,Ϝ,fsl+ϙfusfsl+ϙradT+ϙconvT,E4

where

  • ϙcondT,ϜϙTi+CϙTiC(Ta)
ϙcondT,Ϝ=ϙTi+C(T)+ϙTiC(T)-ϙTi+C(T).Ϝ.fT,E5

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

  • ϙfusfsl=fsl1-ϜT=Tsl=1943fsl1-Ϝ
  • ϙradTdpϓp
ϙradT=4ϓϡT3dpp-3.E6
  • ϙconvTϙairT
ϙconvT=ϙairTp-3.E7

2.2.3. Enthalpy balance

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

ϟdhdt=.ϙT,E8

where ϟin kg.m-3is the density of the system, hthe mass enthalpy and ϙthe thermal conductivity. Temperature Tand 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Ϗϐ,fslthat describe the evolution of the system. The mass enthalpy is defined by

hT,Ϝ,fϏϐ,fsl=ϜβrHTa+TaTCpϖ,Ϝdϖ+fϏϐLϏϐ+fslLsl.E9

The computation of the partial derivatives hT, hϜ, hfϏϐand hfslleads to

hT=CpT,Ϝ,hfϏϐ=LϏϐ,hfsl=Lsl,hϜ=βrHTa+TaTϜCpϖ,Ϝϖ.E10

Thanks to the definition of CpT,Ϝgiven by Eq.(3), we obtain

hϜ=βrHTa+TaTCpTiCϖ-MTi.CpTiϖ+MC.CpCϖMTiCdϖ,E11

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

IT=TaTCpTiCϖ-MTiCpTiϖ+MCpCϖMTiCdϖ.E12

Fig.1 represents the evolution of I(T)βrHfor T[300,3500]. It is observed that max300T3500ITβrHTa2.510-5. The contribution of TaTϜCpϖ,Ϝdϖto hϜis therefore neglected.

Finally the enthalpy balance for the system is written by

ϟ.hsT,Ϝt=.ϙT+ϟ.-βrHTa.Ϝt-ϟ.LϏϐ.fϏϐt-ϟ.Lsl.fslt,E13

where

hsT,Ϝ=TaTCpϖ,Ϝϖ.E14

Radiative boundary conditions are defined over χfor the system and take the form

-ϙTn=ϓ.ϡ.T4-Tχ4,E15

where ϓ=0.7is 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.

TM,0=Ta.E16

Figure 1.

Evolution ofI(T)βrHfor T ∈[300, 3500] K.

2.2.4. Adiabatic temperature of the reaction Tad

The adiabatic temperature Tadof 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 βrHTa=βfH0TiCTa. Considering a thermodynamical path composed of the two steps:

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

  2. End product TiC is heated up from Tato Tad.

One obtains that

βfH0TiCTa+TaTadCpTiCdT=0.E17

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

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.

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 χi×tn,tn+1, where χi=[ri,ri+1]is a control volume. A cell-centered approximation is used. The unknown Tindenotes the mean value of T(x,t)at time tnover χi. mi=riri+1rgdris the discrete "mass" of control volume χi.

tntn+1χiϟ.hsT,Ϝtdχdt=tntn+1χiϙT,Ϝ,fsl+ϟ.βHfMTiCkT1-Ϝdχdt.E18

Backward Euler implicit scheme is used for the temporal discretization since it is unconditionally stable, robust and well adapted for the discretization of such problems. βtis 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

miϟhsTin+1,Ϝin+1-hsTin,Ϝinβt=ϙi+12n+1ri+1gTi+1n+1-Tin+1di+12-ϙi-12n+1rigTin+1-Ti-1n+1di-12+miϟβHfMTiCkd(Tn+1)1-Ϝin+1.E19

Here di+12is the distance between the center of two adjacent control volumes, and denoting by hi=ri+1-rithe length of control volume χi, then di+12=(hi+1+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, ϙi+12n+1is computed thanks to the following harmonic mean formula:

ϙi+12n+1=hi+hi+1hi+1ϙi+1n+1+hiϙin+1,E20

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

3.1.2. Numerical discretization of the exothermic kinetics

We integrate the differential equation over a space-time finite volume χi×tn,tn+1. A cell-centered approximation is used.

tntn+1χidϜdt=tntn+1χikdT1-Ϝ.E21

The backward Euler scheme, first order accurate fully implicit scheme is used. Denoting by Ϝinthe mean value at time tnover χiof Ϝ(x,t), we obtain the following discrete equation

Ϝin+1=Ϝin+βtkd(Tin+1)1+βtkd(Tin+1).E22

It was proved in [1] that the numerical scheme is Lstable, moreover for each time index nand for all control volume index i=1,,I:Ϝin[0,1]. For a given control volume index i, the time sequence (Ϝin)nis 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,j×tn,tn+1. Each rectangular structured non uniform control volume χi,j=[ri,ri+1]×[zj,zj+1], has surface mi,jdefined by :

mi,j=χi,jdχ=zjzj+1riri+1rgdz=ri+1g+1-rig+1g+1.zj+1-zj.E23

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

mi,jβt.ϟ.hsTi,jn+1,Ϝi,jn+1-hsTi,jn,Ϝi,jn=ϙi+12,jn+1ri+1gTi+1,jn+1-Ti,jn+1di+12r-ϙi-12,jn+1rigTi,jn+1-Ti-1,jn+1di-12r+ϙi,j+12n+1Ti,j+1n+1-Ti,jn+1dj+12z-ϙi,j-12n+1Ti,jn+1-Ti,j-1n+1dj-12z+mi,j.ϟ.βHfMTiCkTi,jn+11-Ϝi,jn+1.

The expression used for ϙi+12,jn+1, ϙi-12,jn+1, ϙi,j+12,jn+1, ϙi,j-12n+1is 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,k×tn,tn+1. Each parallelepipedic structured non uniform control volume χi,j,k=[xi,xi+1]×[xj,yj+1]×[zk,zk+1], has volume mi,j,k=χi,j,kdχ.

mi,j,kβt.ϟ.hsTi,j,kn+1,Ϝi,j,kn+1-hsTi,j,kn,Ϝi,j,kn=ϙi+12,j,kn+1Ti+1,j,kn+1-Ti,j,kn+1di+12x-ϙi-12,j,kn+1Ti,j,kn+1-Ti-1,j,kn+1di-12x+ϙi,j+12,kn+1Ti,j+1,kn+1-Ti,j,kn+1dj+12y-ϙi,j-12,kn+1Ti,j,kn+1-Ti,j-1,kn+1dj-12y+ϙi,j,k+12n+1Ti,j,k+1n+1-Ti,j,kn+1dk+12z-ϙi,j,k-12n+1Ti,j,kn+1-Ti,j,k-1n+1dk-12z+mi,j,k.ϟ.βHfMTiCkTi,j,kn+11-Ϝi,j,kn+1.E24

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

3.4. Numerical implementation

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

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

In this section we analyze the ignition and propagation of the combustion front during titanium carbide-TiC- combustion synthesis. It is worth mentioning that the propagation of the combustion wave is either stable i.e. at constant velocity or oscillatory and can be determined according to [16] by the computation of the Zeldovich number Zedefined by Ze=Tad-TaE*Rad2. Tad[K]is the adiabatic temperature of the reaction, Ta[K]is the ambiant temperature and kd*[s-1],E*[kJ/mol]characterizes the first order exothermic kinetics. Ris 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 (k0,E*)=(3800s-1,30.0kcal.mol-1), the Zeldovich number Ze=4.52<Zec=22+5, therefore the solid flame propagation is stable.

4.1. 1D numerical study

This subsection presents unpublished results related to the induction time Ϣindin 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 Reare

-ϙTn0,t=ϓϡT0,t4-Tfx=04,-ϙTnRe,t=ϓϡTRe,t4-Tfx=Re4.E25

4.1.1. Tools used to analyse the numerical simulation results

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

  • The induction time Ϣindis 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 xind>0such that Ϝ(xind,Ϣind)=0.5,

  • Ending time of the reaction is time Ϣend>0such that xχ,Ϝ(x,Ϣend)=1,

  • 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 computation of spatial and temporal evolution of points xM0.01(t), xM0.50(t)et xM0.99(t)such that Ϝ(xM0.01(t),t)=0.01, Ϝ(xM0.50(t),t)=0.50et Ϝ(xM0.99(t),t)=0.99,

  • A each instant t, xmax(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

Tsyn(x)=0+T(x,t)dϜdt(x,t)dx.E26

It was defined and used in [2].

A meaningful numerical simulation is presented in Fig.(2) and in Fig.(3). We observe that the "thickness" of the reaction-zone defined by xM0.99(t)-xM0.01(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=0=1600K. 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 xM0.01(t),xM0.99(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.

Figure 2.

Temporal evolution of xmax(t) and x(Ϝ=0.01,t) x(Ϝ=0.50,t) x(Ϝ=0.99,t) for wall temperature Tfr=0=1600 K. Phase changes are not taken into account.

Figure 3.

Temporal evolution of xmax(t) and x(Ϝ=0.01,t) x(Ϝ=0.50,t) x(Ϝ=0.99,t) for wall temperature Tfr=0=1600 K. Phase changes are taken into account.

4.1.2. Combustion front propagation

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

Figure 4.

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

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

Figure 5.

Temporal evolution of ϟβHfMTiCkd(T)1-Ϝi at equally spaced instants. Phase change is not taken into account.

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 Ϣindand ending time Ϣendincrease exponentially when Tfdecreases, whenever phase change is taken into account or not.

Figure 6.

Evolution of Ϣind et Ϣend as a function of temperature furnace Tfr=0.

Fig.(7) shows that the evolution of xind(t)is nearly similar whatever the furnace temperature Tfis. 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.

Figure 7.

Temporal evolution of xind(t) position, for various values of the furnace temperature Tfr=0∇800,1200,1600,2000,2400K.

Figure 8.

Time evolution of Tmax(t) for several values of temperature furnace Tfr=0∇800,1200,1600,2000,2400.

Fig.(9) shows the same behaviour for xmax(t)whatever Tfis. 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 xmax(t)=0, so the temperature is below the ignition temperature.

  • A combustion front propagation step xmax(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=Re. 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 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=0doesn’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.

Figure 9.

Time evolution of xmax(t) for different values of furnace temperature Tfr=0.

Spatial distribution for different values of Tfr=0for 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.

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.

4.1.4. Energy released/absorbed by the boundary conditions

We analyze the time evolution of flux τ0t, exchanged at r=0between the exterior and the material, defined by τ0t=-ϙTn0,t=ϓϡT(0,t)4-Tf4. Three consecutive steps can be observed thanks to the analysis of τ0trepresented by Fig.(12) and correlated to the time-evolution of T(0,t)

Figure 12.

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

  1. T(0,t)TaTfτ0t
  2. T(0,t)TfTadTf<Tadτ0t
  3. T(0,t)Tfτ0t

A similar analysis can be done for τRet=-ϙTnRe,t=ϓϡT(Re,t)4-Tf4. It is worth mentioning that due to the thermal shock observed when the "hot" solid combustion front reaches the "cold" boundary, the amplitude τRetis an order of magnitude higher than τ0t. 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 Tsequal to the first phase transition TiϏTiϐtemperature, TϏϐ=1166K. We observe on Fig.(13), that the discrepancy ϢignTs=0-ϢignTs=TϏϐincreases significantly when Tfdecreases whenever the phase change is taken into account or not. In practice, ϢignTs=0Tf=900K=ϢignTs=TϏϐTf=1600K. A similar conclusion is drawn on Fig.(14) for the difference ϢendTs=0-ϢendTs=TϏϐ.

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

Figure 13.

Evolution of Ϣind as a function of tfr=0 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=0 for a kinetics with/without a cut-off temperature. Phase change taken into account (PC=1), or not (PC=0).

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

Figure 15.

Spatial temperature profile at t=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

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

Figure 17.

Spatial temperature profiles at several consecutive instants.

Figure 18.

Spatial conversion rate profiles at several consecutive instants.

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

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.

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 Ϣindslab>Ϣindcyl>Ϣindsph, for each value of Tf. Curvature effects may explain this result.

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

We analyze the contribution of the exponent ddefined in Eq.(2) to the induction time Ϣindand 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 k0(T)=k0*exp-E*RTand kd(T)=kd*Tdexp-E*RTat T=Tad

k0*=kd*.Tadd.E27

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 dwe take into account or not the phase change and use Tfr=0=1600K.

Figure 22.

Temperature dependance of kd(T) for various values of exponent d∇[0,2].

We point out that

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

  • When TTad, then kd(T)d>0is below k0(T),

  • When T>Tad, then kd(T)>>k0(T). This imply that the heat released by the exothermic kinetics is significantly increasing when dincreases. Temperature "overshoot" of Tmax(t)can be observed on Fig.(30).

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

Fig.(23) shows that the increasing rate of the induction time Ϣindis super-linear, while it is linear for the ending time Ϣendwhen degree dincreases. It appears independent from the eventual phase change contribution.

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=0=1600 K and Tfr=Re=300 K.

It is worth mentioning that all spatial temperature profiles represented on Fig.(24) have a discontinuity of the time-derivative when T=Tsl. This is explained by the contribution ϙfusfslof 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=0=1600 K and Tfr=Re=300 K.

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

Figure 25.

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

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

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=0=1600 K and Tfr=Re=300 K.

Fig.(27) shows that the time evolution of xmax(t)has a similar shape for each value of dphase change is taken into account or not. Considering high values of dimply 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 Ϣindwhen high values of dare 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 dsince 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 das seen on Fig.(30) is in fact nearly independent from dwhen a suitable time-translation is performed whenever phase change is taken into account or not.

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=0=1600 K and Tfr=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=0=1600 K and Tfr=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=0=1600 K and Tfr=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=0=1600 K and Tfr=Re=300 K.

4.2. 2D numerical study

This subsection computes the propagation of a dual radial/longitudinal front and a single longitudinal front for various values of the temperature furnace in which the cylindrical sample is placed. The main interest of this modelling with respect to the 1D case is to analyze the contribution of the lateral heat losses over the ignition and propagation of the combustion front inside the cylindrical sample. The boundary conditions are defined over χ=αr=0αr=Reαz=0αz=Heby

r,zαr=0:-ϙT0,z,tnr=0e=0,r,zαr=Re:-ϙTRe,z,tnr=Ree=ϓ.ϡ.TRe,z,t4-Tfz,r=Re(t)4,r,zαz=0:-ϙTr,0,tnz=0e=ϓ.ϡ.Tr,0,t4-Tfr,z=0t4,r,zαz=He:-ϙTr,He,tnz=Hee=ϓ.ϡ.Tr,He,t4-Tfr,z=He(t)4.E28

The choice of the temperature triple Tfz,r=Re(t),Tfr,z=0(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=300K. Each computation will analyze the phenomena during t=10s. The computational grid is uniform along z-axis and iso-volume along r-axis.

4.2.1. Reaction-diffusion vs thermal diffusion

We choose Tfz,r=Re(t),Tfr,z=0(t),Tfr,z=He(t)=300K,1600K,300Kfor 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 TaT(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.

Figure 31.

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

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

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=0(t),Tfr,z=He(t)=2400K,2400K,300K. 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 32.

Temperature distribution T(r,z,3) for (r,z)∇χ.

Figure 33.

Integral of the energy released by the kinetics during the process at each (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 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)∇χ.

Figure 35.

Time distribution of the thermal history time T1800K(r,z,5) for (r,z)∇χ.

Figure 36.

Time distribution of the thermal history time T2400K(r,z,5) for (r,z)∇χ.

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

Figure 37.

Temperature distribution T(x,y,z,3) for (x,y,z)∇χ.

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.

© 2012 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

A. Aoufi and G. Damamme (September 19th 2012). Multidimensional Numerical Simulation of Ignition and Propagation of TiC Combustion Synthesis, Numerical Simulation - From Theory to Industry, Mykhaylo Andriychuk, IntechOpen, DOI: 10.5772/48441. Available from:

chapter statistics

1085total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Numerical Simulation of Combustion in Porous Media

By Arash Mohammadi and Ali Jazayeri

Related Book

First chapter

Modeling the Physical Phenomena Involved by Laser Beam – Substance Interaction

By Marian Pearsica, Stefan Nedelcu, Cristian-George Constantinescu, Constantin Strimbu, Marius Benta and Catalin Mihai

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More about us