adiabatic Simulation Using TTv model

[sramjatan1]

Hi, I am trying to run a simple adiabatic simulation of O2 dissociation using two temperatures,
starting with the example in /c++/O2_dissociation/O2_dissociation.cpp.

The example file which comes with the mutation distribution simulates a constant volume, adiabatic reactor beginning
with pure O-atoms at 4000 K and a density of 2.85e-4 kg/m^3. However, the ChemNonEq1T model is used.
I would like to use the ChemNonEqTTv model instead. My c++ is not too good. I would appreciate any insight on getting this example to work. I have a feeling I am not calling the setState routine correctly. I tried to make a vector T_test which holds the initial translational (2000 K), and vibrational (300 K) temperatures.

The code is directly taken from the examples folder and hence I am not including the printHeader and printResults functions in the code below.

Example file

Code

**Code**
```c++


// Main entry point
int main()
{
    // Initial conditions are defined here
    const double T_init   = 4000.0;  // K
    const double rho_init = 2.85e-4; // kg/m^3

    vector<double> T_test(20000.0,1000.0);
 

    // First create the mixture object
    MixtureOptions opts;
    opts.setSpeciesDescriptor("O O2");        // 2 species O2 mixture
    opts.setThermodynamicDatabase("RRHO");    // Thermo database is RRHO
    opts.setMechanism("O2");                  // O2 + M = 2O + M
  //  opts.setStateModel("ChemNonEq1T");        // chemical noneq. w/ 1T
    opts.setStateModel("ChemNonEqTTv");        // chemical noneq. w/ 2T


    Mixture mix(opts);                        // Init. the mixture with opts

    // Setup arrays
    double rhoi [mix.nSpecies()];
    double wdot [mix.nSpecies()];

    // Set state of mixture to the initial conditions
    rhoi[mix.speciesIndex("O")]  = 0.0;
    rhoi[mix.speciesIndex("O2")] = rho_init;
   // mix.setState(rhoi, &T_init, 1); // set state using {rho_i, T}

    mix.setState(rhoi, T_test.data(), 1);
 

    // Get the mixture energy which is conserved during the simulation
    double etot = mix.mixtureEnergyMass() * mix.density(); // J/m^3

    // Write the results header and initial conditions
    printHeader(mix);
    printResults(0.0, mix);

    // Integrate the species mass conservation equations forward in time
    double dt   = 0.0;
    double time = 0.0;
    double tout = 0.0001;

    while (time < 100.0) {
        // Get the species production rates
        mix.netProductionRates(wdot);

        // Compute "stable" time-step based on maximum allowed change in
        // species densities
        dt = 0.0;
        for (int i = 0; i < mix.nSpecies(); ++i)
            dt += 1.0e-7 * std::abs(rho_init * mix.Y()[i] / wdot[i]);
        dt = min(dt / mix.nSpecies(), tout-time);

        // Integrate in time
        for (int i = 0; i < mix.nSpecies(); ++i)
            rhoi[i] += wdot[i]*dt;
        time += dt;

        // Set the new state of the mixture (using conserved variables)
        mix.setState(rhoi, &etot);

        // Print results at specified time intervals
        if (std::abs(tout-time) < 1.0e-10) {
           printResults(time, mix);

           
           if (tout > 9.9)          tout += 10.0;
           else if (tout > 0.99)    tout += 1.0;
           else if (tout > 0.099)   tout += 0.1;
           else if (tout > 0.0099)  tout += 0.01;
           else if (tout > 0.00099) tout += 0.001;
           else                     tout += 0.0001;
       
           }
    }

    // Now get the equilibrium values
    mix.equilibriumComposition(mix.T(), mix.P(), rhoi); // puts X_i in rhoi
    mix.convert<X_TO_Y>(rhoi, rhoi);                    // converts X_i to Y_i

    cout << "Equilibrium mass fractions at " << mix.T() << " K and "
         << mix.P() << " Pa:" << endl;

    for (int i = 0; i < mix.nSpecies(); ++i)
        cout << setw(15) << mix.speciesName(i) << ":"
             << setw(15) << rhoi[i] << endl;

    return 0;
}

Comments
If there are any working examples for a two-temperature adiabatic simulation, please let me know. Thanks!!