Catalytic combustion of a stagnation flow on a platinum surface

This script solves a catalytic combustion problem. A stagnation flow is set up, with a gas inlet 10 cm from a platinum surface at 900 K. The lean, premixed methane/air mixture enters at ~ 6 cm/s (0.06 kg/m2/s), and burns catalytically on the platinum surface. Gas-phase chemistry is included too, and has some effect very near the surface.

The catalytic combustion mechanism is from Deutschmann et al., 26th Symp. (Intl.) on Combustion,1996 pp. 1747-1754

Tags: Matlab combustion catalysis 1D flow surface chemistry

Initialization

help catcomb;

clear all
close all

t0 = cputime; % record the starting time

Set parameter values

p = OneAtm; % pressure
tinlet = 300.0; % inlet temperature
tsurf = 900.0; % surface temperature
mdot = 0.06; % kg/m^2/s
transport = 'mixture-averaged'; % transport model

Solve first for a hydrogen/air case for use as the initial estimate for the methane/air case.

% composition of the inlet premixed gas for the hydrogen/air case
comp1 = 'H2:0.05, O2:0.21, N2:0.78, AR:0.01';

% composition of the inlet premixed gas for the methane/air case
comp2 = 'CH4:0.095, O2:0.21, N2:0.78, AR:0.01';

% the initial grid, in meters. The inlet/surface separation is 10 cm.
initial_grid = [0.0, 0.02, 0.04, 0.06, 0.08, 0.1]; % m

% numerical parameters
tol_ss = {1.0e-8 1.0e-14}; % {rtol atol} for steady-state problem
tol_ts = {1.0e-4 1.0e-9}; % {rtol atol} for time stepping

loglevel = 1; % amount of diagnostic output (0 to 5)

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

Create the gas object

This object will be used to evaluate all thermodynamic, kinetic, and transport properties

The gas phase will be taken from the definition of phase gas in input file ptcombust.yaml, which is a stripped-down version of GRI-Mech 3.0.

gas = Solution('ptcombust.yaml', 'gas', transport);
gas.TPX = {tinlet, p, comp1};

Create the interface object

This object will be used to evaluate all surface chemical production rates. It will be created from the interface definition Pt_surf in input file ptcombust.yaml, which implements the reaction mechanism of Deutschmann et al., 1995 for catalytic combustion on platinum.

surf_phase = Interface('ptcombust.yaml', 'Pt_surf', gas);
surf_phase.TP = {tsurf, surf_phase.P};

Integrate the coverage equations in time for 1 s, holding the gas composition fixed to generate a good starting estimate for the coverages.

surf_phase.advanceCoverages(1.0);

The two objects we just created are independent of the problem type – they are useful in zero-D simulations, 1-D simulations, etc. Now we turn to creating the objects that are specifically for 1-D simulations. These will be ‘stacked’ together to create the complete simulation.

Create the flow object

The flow object is responsible for evaluating the 1D governing equations for the flow. We will initialize it with the gas object, and assign it the name flow.

flow = AxisymmetricFlow(gas, 'flow');

% set some parameters for the flow
flow.P = p;
flow.setupGrid(initial_grid);
flow.setSteadyTolerances('default', tol_ss{:});
flow.setTransientTolerances('default', tol_ts{:});

Create the inlet

The temperature, mass flux, and composition (relative molar) may be specified. This object provides the inlet boundary conditions for the flow equations.

inlt = Inlet(gas, 'inlet');

% set the inlet parameters. Start with comp1 (hydrogen/air)
inlt.T = tinlet;
inlt.massFlux = mdot;
inlt.setMoleFractions(comp1);

Create the surface

This object provides the surface boundary conditions for the flow equations. By supplying object surface_phase as an argument, the coverage equations for its surface species will be added to the equation set, and used to compute the surface production rates of the gas-phase species.

surf = ReactingSurface(surf_phase, 'surface');
surf.T = tsurf;

Create the stack

Once the component parts have been created, they can be assembled to create the 1D simulation.

stack = Sim1D({inlt, flow, surf});

% set the initial profiles.
stack.setProfile(2, {'velocity', 'spread_rate', 'T'}, ...
                [0.0, 1.0 % z/zmax
                 0.06, 0.0 % u
                 0.0, 0.0 % V
                 tinlet, tsurf]); % T
names = gas.speciesNames;

for k = 1:gas.nSpecies
    y = inlt.massFraction(k);
    stack.setProfile(2, names{k}, [0, 1; y, y]);
end

stack.setTimeStep(1.0e-5, [1, 3, 6, 12]);
stack.setMaxJacAge(4, 5);

Solution

Start with the energy equation on

flow.energyEnabled = true;

Disable the surface coverage equations, and turn off all gas and surface chemistry

surf.coverageEnabled = false;
surf_phase.setMultiplier(0.0);
gas.setMultiplier(0.0);

Solve the problem, refining the grid if needed

stack.solve(1, refine_grid);

Now turn on the surface coverage equations, and turn the chemistry on slowly

surf.coverageEnabled = true;

for iter = 1:6
    mult = 10.0^(iter - 6);
    surf_phase.setMultiplier(mult);
    gas.setMultiplier(mult);
    stack.solve(1, refine_grid);
end

At this point, we should have the solution for the hydrogen/air problem. Now switch the inlet to the methane/air composition.

inlt.setMoleFractions(comp2);

Set more stringent grid refinement criteria

stack.setRefineCriteria(2, 100.0, 0.15, 0.2);

Solve the problem for the final time

stack.solve(loglevel, refine_grid);

Show statistics

stack.writeStats;
elapsed = cputime - t0;
e = sprintf('Elapsed CPU time: %10.4g', elapsed);
disp(e);

Make plots

clf;

subplot(3, 3, 1);
plotSolution(stack, 'flow', 'T');
title('Temperature [K]');

subplot(3, 3, 2);
plotSolution(stack, 'flow', 'velocity');
title('Axial Velocity [m/s]');

subplot(3, 3, 3);
plotSolution(stack, 'flow', 'spread_rate');
title('Radial Velocity / Radius [1/s]');

subplot(3, 3, 4);
plotSolution(stack, 'flow', 'CH4');
title('CH4 Mass Fraction');

subplot(3, 3, 5);
plotSolution(stack, 'flow', 'O2');
title('O2 Mass Fraction');

subplot(3, 3, 6);
plotSolution(stack, 'flow', 'CO');
title('CO Mass Fraction');

subplot(3, 3, 7);
plotSolution(stack, 'flow', 'CO2');
title('CO2 Mass Fraction');

subplot(3, 3, 8);
plotSolution(stack, 'flow', 'H2O');
title('H2O Mass Fraction');

subplot(3, 3, 9);
plotSolution(stack, 'flow', 'H2');
title('H2 Mass Fraction');