sandbox/ecipriano/src/reactors.h

Batch Reactors

Collection of ideal reactor functions for the operator splitting solution of reactive systems (Cuoci et al. 2013). The temperature and chemical species concentrations are considered perfectly mixed inside each cell. Neglecting the convective and diffusive transport in the conservation equations, the equations reduce to the ODE system of a batch reactor: \displaystyle \begin{cases} \dfrac{d\omega_i}{dt} = \dfrac{1}{\rho} \sum \limits_{j=1}^{NR} \nu_{i,j} r_j MW_i \\ \dfrac{d T}{dt} = \dfrac{1}{\rho c_p} \sum \limits_{j=1}^{NR} r_j \Delta H_{R_j} \end{cases} where \nu is the matrix of the stoichiometric coefficients, while r_j is the rate of the j-th reaction (Fogler 2020). The term \Delta H_R is the enthalpy change associated with the chemical reactions. The reaction rate r_j is usually computed as: \displaystyle r_j = k_j \prod \limits_{i=1}^{NS} c_i^{\nu_{i,j}} and the kinetic constant k_j is computed from a modified Arrhenius equation: \displaystyle k_j = AT^n \exp^{-E_a/(RT)} This procedure is managed by the OpenSMOKE++ library (opensmoke.h).

#include "radiation.h"

User Data

Structure that gathers arguments that can be used inside the batch reactor system functions.

typedef struct {
  double rho;
  double cp;
  double P;
  double T;
} UserDataODE;

batch_isothermal_constantpressure(): solve chemical species reactions

• y: vector (length = NS) containing the initial values of the system
• dt: simulation time step
• dy: right-hand-side of the batch system of equations
• args: structure with additional arguments to be passed to the system

void batch_isothermal_constantpressure (const double * y, const double dt, double * dy, void * args) {

Unpack data for the ODE system.

  UserDataODE data = *(UserDataODE *)args;
  double rho = data.rho;
  double cp = data.cp;
  double Pressure = data.P;
  double Temperature = data.T;

Set temperature and pressure of the system, for the calculation of the kinetic constants.

  OpenSMOKE_GasProp_SetTemperature (Temperature);
  OpenSMOKE_GasProp_SetPressure (Pressure);

We create a vector with the mass fractions and we convert it to mole fractions.

  double massfracs[NGS], molefracs[NGS];
  for (int jj=0; jj<NGS; jj++) {
    massfracs[jj] = y[jj];
  }

  mass2molefrac (molefracs, massfracs, inMW, NGS);
  double MWMix = mass2mw (massfracs, inMW, NGS);

  double ctot = Pressure/(R_GAS*1000.)/Temperature;
  double ci[NGS], ri[NGS];
  for (int jj=0; jj<NGS; jj++) {
    ci[jj] = ctot*molefracs[jj];
    ri[jj] = 0.;
  }
  rho = ctot*MWMix;

Compute the kinetic constant, reaction rates, and formation rates.

  OpenSMOKE_GasProp_KineticConstants ();
  OpenSMOKE_GasProp_ReactionRates (ci);
  OpenSMOKE_GasProp_FormationRates (ri);

Equation for the chemical species.

  for (int jj=0; jj<OpenSMOKE_NumberOfSpecies(); jj++) {
    dy[jj] = OpenSMOKE_MW(jj)*ri[jj]/rho;
  }
}

batch_nonisothermal_constantpressure(): solve chemical species reactions

• y: vector (length = NS+1) containing the initial values of the system
• dt: simulation time step
• dy: right-hand-side of the batch system of equations
• args: structure with additional arguments to be passed to the system

void batch_nonisothermal_constantpressure (const double * y, const double dt, double * dy, void * args) {

Unpack data for the ODE system.

  UserDataODE data = *(UserDataODE *)args;
  double rho = data.rho;
  double cp = data.cp;
  double Pressure = data.P;
  double Temperature = y[OpenSMOKE_NumberOfSpecies()];

Set temperature and pressure of the system, for the calculation of the kinetic constants.

  OpenSMOKE_GasProp_SetTemperature (Temperature);
  OpenSMOKE_GasProp_SetPressure (Pressure);

We create a vector with the mass fractions and we convert it to mole fractions.

  double massfracs[NGS], molefracs[NGS];
  for (int jj=0; jj<NGS; jj++) {
    massfracs[jj] = y[jj];
  }

  mass2molefrac (molefracs, massfracs, inMW, NGS);
  double MWMix = mass2mw (massfracs, inMW, NGS);

  double ctot = Pressure/(R_GAS*1000.)/Temperature;
  double ci[NGS], ri[NGS];
  for (int jj=0; jj<NGS; jj++) {
    ci[jj] = ctot*molefracs[jj];
    ri[jj] = 0.;
  }
  rho = ctot*MWMix;

  OpticallyThinProperties otp;
  otp.T = Temperature;
  otp.P = Pressure;
  otp.xH2O = molefracs[otm.indexH2O];
  otp.xCO2 = molefracs[otm.indexCO2];

Compute the kinetic constant, reaction rates, and formation rates.

  OpenSMOKE_GasProp_KineticConstants ();
  OpenSMOKE_GasProp_ReactionRates (ci);
  OpenSMOKE_GasProp_FormationRates (ri);

Equation for the chemical species.

  for (int jj=0; jj<OpenSMOKE_NumberOfSpecies(); jj++) {
    dy[jj] = OpenSMOKE_MW(jj)*ri[jj]/rho;
  }

Get the heat of reaction and compute the equation for the temperature. We add non-linear contributions such as the heat dissipated for radiation

  double QR = OpenSMOKE_GasProp_HeatRelease (ri);
  dy[OpenSMOKE_NumberOfSpecies()] = (QR + divq_rad (&otp))/rho/cp;
}

References