OpenMP parallel computing of 2D TiC combustion synthesis process using an explicit finite-volume scheme

This paper analyzes from a numerical point of view the ignition and propagation of the combustion front during the exothermic TiC combustion synthesis of a material made of pressed titanium and carbide particles when thermophysical properties are either assumed constant or temperature and conversion rate dependent. A two-dimensional cylindrical geometry is considered. The heat supply is prescribed on one, two or three sides of the physical domain. A one-step kinetics is used to describe the reaction Ti+C→TiC in a solid phase and leads to the computation of the conversion rate. A coupling with a non-linear heat equation which takes into account the heat generated by the exothermic kinetics and the two allotropic phase-changes is considered. An explicit finite-volume discretization of the coupled system is constructed and analyzed. Time-step’s stability condition is given for a general expression of the thermo-physical characteristics. A discrete maximum principle is reported. Open MP API was used to parallelize the numerical software written in C. An average speedup of three was obtained on an intel quad-core processor i7-2600. The ignition time and the fraction of unreacted material are systematically computed and compared for several heat supply scenario.


Introduction
The present paper is devoted to the numerical computation of TiC combustion synthesis [1,2] in a 2D cylindrical coordinate system using an explicit scheme implemented on a multicore architecture using OpenMP API. It is organized as follows. Section two reports the main features of the mathematical modelling of this combustion synthesis process (coupling a partial differential equation with a differential equation) and the expression used for thermophysical parameters (thermal conductivity and heat capacity). Section three outlines the explicit finite-volume discretization of both equations and reports a discrete maximum principle for both equations. Section four explains the main ingredients of the parallelization strategy based upon inserting OpenMP #pragma directives and analyzes the efficiency of the software on a multicore intel i7-2600 cpu. Section five presents some numerical simulation results used to systematically investigate the evolution of ignition curves when thermophysical properties are either constant or temperature and conversion rate dependent. The stability of the combustion wave is investigated thanks to the Zeldovitch number. The contribution of allotropic phase changes is also taken into account numerically. The effect of heat supply on one, two and three faces of the cylinder is also reported. A conclusion summarizes the results presented in this study.

Mathematical modeling
Temperature T in K, and conversion rate ∈ [0,1] are the key variables used to describe the evolution of the exothermic equation Ti+C →TiC. The mathematical modeling couples the enthalpy balance with a differential equation.

Nonlinear reaction diffusion and boundary condition
An enthalpy balance written for a cylinder Ω with radius R and height H, leads to the following nonlinear reaction-diffusion equation (1) Along the symmetry axis of the cylinder, the adiabatic boundary condition is prescribed (2) Radiative boundary conditions are defined for the three remaining subparts Γ of ∂Ω. The value of T∞,Γ ranges from 300 K ( room temperature) to 2400 K. is the emissivity of the material and is the Stefan-Boltzmann constant.

Heat capacity
Assuming that the mass of the system remains constant, a linear mixing law between the heat capacity of reactants CpTi(T), CpC(T) [3] and the heat capacity of the product CpTiC (T ) [3] weighted by the conversion rate of the reaction gives ( Figure 1)

Analysis of the numerical simulations
The following three quantities are systematically computed in order to investigate the results of the numerical simulations presented in next sections: • the synthesis temperature Tsyn(x) ( ∈ Ω) is the weighted mean of the temperature by the kinetics defined (from [7]) by ∀ ∈ Ω, • the normalized fraction of unreacted material nfum(t) is defined (from [8]) • the maximum of the temperature Tmax(t) and the mean temperature ̅ ( ) of T (., t) over Ω are computed at each time t > 0.

Finite-volume discretization
The computational domain Ω is decomposed into a set of NI× NJ rectangular control volumes Ωi,j with "mass" mi,j . A conservative finite-volume cell-centered approximation is used to discretize the governing equations of the modeling. The handling of each phase change transition during which a liquid fraction appears is done thanks to a specific procedure.

Explicit finite-volume discretization of the reaction-diffusion equation
The integration of equation (1) Under the assumption that 0 < ∆t ≤ ∆tstab, the following discrete maximum principle properties for the time-explicit discretization ( [7], [9] for the time-implicit discretization) are established by induction It is worth mentioning that when the thermo-physical properties are constant, the time-step for the reaction-diffusion is also constant, therefore the stability condition Δ is computed once and passed as a parameter to the function advancing the heat transfer the case. When these parameters are temperature and conversion rate dependent, Δ is recomputed at each step once the positive quantities , , , , , , , are evaluated as a function of T and . The computational complexity is therefore a multiple of × .

Choice of the time-step
The time-step ∆tn is computed thanks to the expression Δt = (Δt , Δt ).
(13) It is worth mentioning that in the numerical simulations presented in next sections, we always observed that Δt < . This is explained by the small value of the activation energy * .   Figure 3 shows the scalability of the numerical code. A speedup of 3 was obtained when the number of threads is 4 for a given mesh. It varies slightly as a function of the mesh. Figure 4 shows that the efficiency of the numerical code is above 0.75 when four cores are used. This means that time spent on transferring datas between cores is much lower than time spent on computations. The implementation is therefore efficient.

Stability analysis of the combustion front
Linear stability analysis (from [10]) is based upon the computation of Zeldovich number Here ≅ 2900 and * ≅ 30 = 125.4 kJ, then = 4.52 < 2(2 + √5) ≅ 8.4723 for which a stable propagation with a constant velocity is obtained. (c-pc-0), b) variables properties without phase change (v-pc-0) or c) variables properties with phase change (v-pc-1) are considered.

Propagation of a longitudinal front
In this case, (Tz=0, Tr=R, Tz=H ) = (Tf , 300, 300), where Tf =800, 1600, 2400. A uniform tensorial grid made of 32 × 512 cells is used for the computations. Results are presented in Figure 5, Figure 6, Figure 7 and show that nonlinear expression for Cp(T, ) and (T, ) have a leading contribution to the celerity of the synthesis.

Propagation of a converging radial front
In this case, (Tz=0, Tr=R, Tz=H) = (300, Tf , 300), where Tf =800 K, 1600 K, 2400 K. A uniform tensorial grid made of 512 × 32 cells is used for the computations. Results are presented in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12. A conclusion similar to the longitudinal case can be established. Phase changes contribution is noticeable graphically.    Figure 17. A similar conclusion can be established. Moreover, along z = H/2 a melting might occur because Tmax(t) is locally greater than 3000 K. Latent heat of phase change is therefore needed to accurately capture the effects of this melting.

Propagation of three converging fronts
In this case, (Tz=0, Tr=R, Tz=H) = (Tf, Tf, Tf), where Tf=800, 1600, 2400. This corresponds to the practical case where the cylindrical sample is placed inside a furnace with uniform temperature Tf . A uniform tensorial grid made of 512 × 512 cells is used for the computations.
Results are presented in Figure 18, Figure 19, Figure 20, Figure 21, Figure 22 (notice the symmetries induced by the symmetries of the boundary conditions). The conclusion is similar to those previously obtained.

Conclusions and perspectives
In the present paper we have shown through numerical simulations (thanks to a tuned parallel software using OpenMP and implementing an explicit finite-volume scheme for which a discrete maximum principle was established) that independently of the way external heat supply is applied on part of the boundary ∂Ω, the expressions used for Cp(T, ) and (T ) have a significant influence on TiC combustion synthesis process, quantified by nfum(t) and Tsyn(x). The contribution of the two allotropic phase changes ranges from negligible to minimal. A similar approach is currently investigated in the case of a cubic sample where the computations are genuinely three dimensional.