Adiabatic, constant pressure reactor

This example illustrates how to use class Reactor for zero-dimensional kinetics simulations. Here the parameters are set so that the reactor is adiabatic and very close to constant pressure.

Requires: cantera >= 3.2.0

Tags: Matlab combustion reactor network ignition delay plotting

function reactor1(g)
    tic
    help reactor1

    if nargin == 1
        gas = g;
    else
        gas = ct.Solution('gri30.yaml', 'gri30', 'none');
    end

    P = 101325.0;

Set the initial conditions

    gas.TP = {1001.0, P};
    nsp = gas.nSpecies;
    xx = zeros(nsp, 1);
    xx(1) = 0.285;
    xx(4) = 0.142;
    xx(48) = 0.573;
    gas.X = xx;

Create a reactor, and insert the gas

    r = ct.zeroD.IdealGasReactor(gas);

Create a reservoir to represent the environment

    a = ct.Solution('air.yaml', 'air', 'none');
    a.TP = {a.T, P};
    env = ct.zeroD.Reservoir(a);

Define a wall between the reactor and the environment and make it flexible, so that the pressure in the reactor is held at the environment pressure.

    w = ct.zeroD.Wall(r, env);

Set expansion parameter. dV/dt = KA(P_1 - P_2)

    w.expansionRateCoeff = 1.0e6;

Set wall area

    w.area = 1.0;

Create a reactor network and insert the reactor:

    network = ct.zeroD.ReactorNet({r});

    nSteps = 100;
    tim(nSteps) = 0;
    temp(nSteps) = 0;
    x(nSteps, 3) = 0;
    t = 0.0;
    dt = 1.0e-5;
    t0 = cputime;

    for n = 1:nSteps
        t = t + dt;
        network.advance(t);
        tim(n) = network.time;
        temp(n) = r.T;
        x(n, 1:3) = r.phase.moleFraction({'OH', 'H', 'H2'});
    end

    disp(['CPU time = ' num2str(cputime - t0)]);

Plot results

    clf;
    subplot(2, 2, 1);
    plot(tim, temp);
    xlabel('Time (s)');
    ylabel('Temperature (K)');
    subplot(2, 2, 2)
    plot(tim, x(:, 1));
    xlabel('Time (s)');
    ylabel('OH Mole Fraction (K)');
    subplot(2, 2, 3)
    plot(tim, x(:, 2));
    xlabel('Time (s)');
    ylabel('H Mole Fraction (K)');
    subplot(2, 2, 4)
    plot(tim, x(:, 3));
    xlabel('Time (s)');
    ylabel('H2 Mole Fraction (K)');

    toc
end