sandbox/ecipriano/src/chemistry.h

    Chemical Reactions

    The chemical reactions are solved using the OpenSMOKE++ libraries and an ODE system of equations. Despite the general formulation of the chemistry, so far it is solved just in the gas phase, focusing on the combustion kinetics. Different types of reactions may require the implementation of new reactor systems (reactors.h).

    #include "opensmoke.h"
    #include "reactors.h"

    We initialize the OpenSMOKE++ ODE solver, which is the solver that we are going to use for the chemistry integration. This solver is specifically conceived for stiff reactive systems, and it will automatically sub-step during the whole time step of the simulation.

    event init (i = 0) {
      OpenSMOKE_InitODESolver ();
    }

    We deallocate the OpenSMOKE++ ODE solver.

    event cleanup (t = end) {
      OpenSMOKE_CleanODESolver ();
    }

    Chemistry Event

    This event assigns the batch reactor function, and it calls the ODE system integrator advancing the chemical species mass fraction and temperature fields according to the presence of chemical reactions.

    Set reactor function and number of equations accordingly.

    #ifdef SOLVE_TEMPERATURE
      odefunction batch = &batch_nonisothermal_constantpressure;
      unsigned int NEQ = OpenSMOKE_NumberOfSpecies() + 1;
    #else
      odefunction batch = &batch_isothermal_constantpressure;
      unsigned int NEQ = OpenSMOKE_NumberOfSpecies();
    #endif
    
      foreach() {
        if (f[] < F_ERR) {

    Gather the initial conditions for the ODE solver.

    #ifdef SOLVE_TEMPERATURE
          double y0ode[NGS+1];
    #else
          double y0ode[NGS];
    #endif
          foreach_elem (YGList, jj) {
            scalar YG = YGList[jj];
            y0ode[jj] = YG[];
          }
    #ifdef SOLVE_TEMPERATURE
          y0ode[NGS] = TG[];
    #endif

    Set additional data to be passed to the ODE system.

          UserDataODE data;
          data.P = Pref;
          data.T = 1200.;
    #ifdef SOLVE_TEMPERATURE
          data.rho = rho2;
          data.cp = cp2;
    #endif
          double sources[NEQ];
          data.sources = sources;

    Solve the ODE system using the OpenSMOKE native ODE solver, specifically conceived for stiff reactive systems.

          OpenSMOKE_ODESolver (batch, NEQ, dt, y0ode, &data);

    Recover the results of the ODE system.

          foreach_elem (YGList, jj) {
            scalar YG = YGList[jj];
            YG[] = y0ode[jj];
    
    #ifdef VARPROP
            scalar DYDt2jj = DYDt2[jj];
            DYDt2jj[] += sources[jj]*cm[];
    #endif
          }
    #ifdef SOLVE_TEMPERATURE
          TG[] = y0ode[NGS];
    # ifdef VARPROP
          DTDt2[] += sources[NGS]*cm[];
    # endif
    #endif
        }
      }
    }