Modelling HMX-Aluminium heterogeneous combustion in Cantera

[abhik....@iitgn.ac.in]

Hi,

I am trying to model HMX-Aluminium heterogeneous combustion using the Flame1 function for burner stabilized flat flame in Matlab. I have created an input file “HMX-Al.xml” which contains all the species thermodynamic and transport properties as well as reaction mechanisms. On running the code I am getting the following error:

***********************************************************************

CanteraError thrown by OneDim::timeStep:

Time integration failed.

***********************************************************************

Error in stack_methods (line 13)

    v = ctmethods(90, n, job, a, b);

Error in Stack/solve (line 14)

stack_methods(s.stack_id, 104, loglevel, refine_grid);

Error in gas_phase_flame_test1 (line 104)

solve(fl, loglevel, refine_grid);

***********************************************************************

Since I am new to Cantera I am unable to understand why it is throwing this error. Currently, I am using Cantera 2.3 version.  In this problem, I am trying to model the combustion of solid aluminium nano particles in a gaseous environment created by the decomposition of HMX. Thus, it should be a heterogeneous reaction. But I could not find a way to this in Matlab (using Cantera) and therefore, in this case I have considered aluminium to be a gas. I have added the required reactions and species details in the input file considering them as gas. My aim is to model aluminium as a solid which is burning and releasing energy in the gaseous environment. The oxidizers present in the gaseous environment are O2, H2O and CO2.  

·         Kindly tell me whether this can be done in Cantera using Matlab and if yes how to do it.

·         Also, how to resolve the above mentioned error.

I have attached the input file "HMX-Al.xml" for reference. Following is the code I am using:

W = [30,44,28,46,296,27]; %Molecular weight of species

r = 0.1E-2; %burning rate

patm = 1;

massfraction_gasinput = [0.05,0.07,0.03,0.05,0.55,0.25];

t0 = cputime;  % record the starting time

 

% parameter values

p          =   patm*oneatm;           % pressure

 

tin        =   680;                   % inlet temperature

 

mdot_o     =   r*1.9*1000;            %  kg/m^2/s

 

rxnmech    =  'HMX-Al.xml';           % reaction mechanism file

 

transport  =  'HMX_mix';              % transport model

 

mass_frac  =  massfraction_gasinput(1:6); % mass fractions of the 5 species

 

comp1      = sprintf('CH2O:%d,N2O:%d,H2CN:%d,NO2:%d,HMX:%d,AL(cr):%d',mass_frac); % gas composition

 

X(1)     =  (massfraction_gasinput(1)/W(1))/(sum(massfraction_gasinput./W)); %CH2O

 

X(2)     =  (massfraction_gasinput(2)/W(2))/(sum(massfraction_gasinput./W)); %N2O

 

X(3)     =  (massfraction_gasinput(3)/W(3))/(sum(massfraction_gasinput./W)); %H2CN

 

X(4)     =  (massfraction_gasinput(4)/W(4))/(sum(massfraction_gasinput./W)); %NO2

 

X(5)     =  (massfraction_gasinput(5)/W(5))/(sum(massfraction_gasinput./W)); %HMX

 

X(6)     =  (massfraction_gasinput(6)/W(6))/(sum(massfraction_gasinput./W)); %AL

 

mole_frac = X(1:6);

 

comp      =  sprintf('CH2O:%d,N2O:%d,H2CN:%d,NO2:%d,HMX:%d,AL(cr):%d',mole_frac);  % gas composition

 

initial_grid = (0.03/patm)*[0.0 0.2 0.4 0.6 0.8 1.0];  % m

 

tol_ss    = [1.0e-5 1.0e-13];       % [rtol atol] for steady-state problem

tol_ts    = [1.0e-4 1.0e-13];       % [rtol atol] for time stepping

 

loglevel  = 1;                      % amount of diagnostic output (0

                                    % to 5)

 

refine_grid = 1;                    % 1 to enable refinement, 0 to

                                    % disable

max_jacobian_age = [5, 10];

 

%%%%%%%%%%%%%%%% create the gas object %%%%%%%%%%%%%%%%%%%%%%%%

 

gas = Solution(rxnmech,transport);

 

% set its state as desired

 

set(gas,'T', tin, 'P', p, 'Y', comp1);

 

lambda_gasphase = thermalConductivity(gas);

 

%%%%%%%%%%%%%%%% create the flow object %%%%%%%%%%%%%%%%%%%%%%%

 

f = AxisymmetricFlow(gas,'flow');

set(f, 'P', p, 'grid', initial_grid);

set(f, 'tol', tol_ss, 'tol-time', tol_ts);

 

 

%%%%%%%%%%%%%%% create the burner %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 

burner = Inlet('burner');

set(burner, 'T', tin, 'MassFlux', mdot_o, 'X',comp);

 

%%%%%%%%%%%%%% create the outlet %%%%%%%%%%%%%%%%%%%%%%%%%%%%

 

outlet = Outlet('out');

 

%%%%%%%%%%%%% create the flame object  %%%%%%%%%%%%

 

fl = flame(gas, burner,f,outlet);

 

setMaxJacAge(fl, max_jacobian_age(1),  max_jacobian_age(2));

 

% solve with fixed temperature profile first

 

solve(fl, loglevel, refine_grid);

 

%%%%%%%%%%%% enable the energy equation %%%%%%%%%%%%%%%%%%%%%

 

enableEnergy(f);

 

setRefineCriteria(fl, 2, 200.0, 0.05, 0.1);

 

solve(fl, loglevel, refine_grid);

 

%%%%%%%%%% show statistics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

writeStats(fl);

elapsed = cputime - t0;

e = sprintf('Elapsed CPU time: %10.4g',elapsed);

disp(e);

[Ray Speth]

Hi,

When working with 1D flames, I cannot recommend highly enough (a) upgrading to Cantera 2.4.0 and (b) switching to the Python interface, which has a much better interface to Cantera’s 1D solver. Specifically, when dealing with convergence issues, I would also suggest looking at some of the ideas in the “Flame Debugging“ Jupyter Notebook, which was presented at our latest workshop.

In this particular case, I would be pretty concerned about the treating aluminum particles as gas phase species. For starters, there are problems with the rate constants of the reactions involving Al. For instance, the reaction entry:

  <!-- reaction 0232    -->
  <reaction reversible="no" id="0232">
  <equation>4 AL + 3 O2 =] 2 AL2O3</equation>
      <rateCoeff>
        <Arrhenius>
           <A units="cm3/mol/s">1646676.606</A>
           <b>0</b>
           <E units="cal/mol">17590.8416</E>
        </Arrhenius>
      </rateCoeff>
      <reactants>AL:4.0 O2:3</reactants>
      <products>AL2O3:2.0</products>
    </reaction>

is a global reaction, going by the stoichiometric coefficients, but it is written as if it is an elementary reaction. In the absence of any other information (reaction orders), Cantera must treat this as an elementary reaction, for which the units of the rate constant are not cm^3/mol/s, but cm^18/mol^6/s, and this discrepancy is certainly leading to the rate of this reaction being calculated incorrectly. I’m not sure where this rate constant came from, but adapting it for this purpose will require some care.

Regards,
Ray

[abhik....@iitgn.ac.in]

Dear Ray,

Thanks for the reply. I have actually tried running the same code in Cantera 2.4.0 in Matlab interface also but getting the same error message there. Coming to the global reactions for Aluminum, I have calculated the values of A and E from a aluminum burning time model. I agree that giving a global reaction without the order means cantera would treat it as an elementary reaction and thus there is unit mismatch. Can you tell me how Cantera calculates the reaction rates and then the net production or destruction rates of a particular species? Any documentation which mentions the stepwise calculation done by cantera after reading the inputs of a particular reaction would be helpful. Also, how the units of A are calculated (as you have mentioned in your reply- cm^18/mol^6/s) if the reactions are to be treated as elementary reaction?

Regards,

Abhik

[Ray Speth]

Abhik,

The forward reaction rate R for an elementary reaction is defined as

where is the forward rate constant, is the molar concentration of species k and is the reactant stoichiometric coefficient for species k. Since the reaction rate has units of kmol/m^3/s and the molar concentration has units of kmol/m^3, the rate constant must have units of (m^3/kmol)^(n-1)/s where n is the sum of the reactant stoichiometric coefficients.

Unfortunately, Cantera does not currently check that the units you specify are correct when converting to its native unit system (SI + kmol), so even if the rate constant you’ve specified was correct in the cm + mol based units, it will be off by many orders of magnitude when converted using the incorrect units.

The best reference I can suggest, which covers much of what is implemented in Cantera, is the book “Chemically Reacting Flow” by Kee, Coltrin, Glarborg and Zhu. Reaction rate calculations are covered in Chapter 12 of the second edition.

Finally, I’ll reiterate the suggestion to move to the Python package for all flame calculations, since the capabilities available make it much more feasible to debug problems with solving 1D systems.

Regards,
Ray