It would be nice if you could provide an example for a non-equilibrium calculation with the two temperature model. It might only be me, but I think the two-temperature model is what really distinguishes Mutation++ from other packages like Cantera (which can do chemical non-eq.). In particular, with the existing examples I cannot infer how I should update the state for the two temperature model (the update of the state with density and total energy does not seem to give physical results).
In particular, I consider an example along the lines of https://github.com/mutationpp/Mutationpp/blob/master/examples/c%2B%2B/O2_dissociation/O2_dissociation.cpp
/// [example_code]
// Must include this header file to use the Mutation++ library
#include "mutation++.h"
using namespace Mutation;
using namespace Mutation::Thermodynamics;
using namespace std;
// Prints the header of the output table
void printHeader_vib(const Mixture& mix)
{
cout << setw(8) << "Time[s]";
cout << setw(12) << "T[K]";
cout << setw(12) << "Tv[K]";
cout << setw(12) << "P[Pa]";
cout << setw(12) << "rho[kg/m^3]";
cout << setw(12) << "e[J/kg]";
for (int i = 0; i < mix.nSpecies(); ++i)
cout << setw(12) << "Y_" + mix.speciesName(i);
cout << endl;
}
// Prints the mixture properties at the given time
void printResults_vib(double time, const Mixture& mix)
{
cout << setw(8) << time;
cout << setw(12) << mix.T();
cout << setw(12) << mix.Tv();
cout << setw(12) << mix.P();
cout << setw(12) << mix.density();
cout << setw(12) << mix.mixtureEnergyMass();
for (int i = 0; i < mix.nSpecies(); ++i)
cout << setw(12) << mix.Y()[i];
cout << endl;
}
// Main entry point
int main()
{
// Initial conditions are defined here
const double T_init = 4000.0; // K
const double Tv_init = 200.0; // K
const double rho_init = 3.76e-3; // kg/m^3
// First create the mixture object
MixtureOptions opts("air_5");
//opts.setStateModel("ChemNonEq1T"); // chemical noneq. w/ 1T
opts.setStateModel("ChemNonEqTTv"); // chemical noneq. w/ 2T
opts.setThermodynamicDatabase("RRHO"); // Thermo database is RRHO
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("NO")] = 0.0;
rhoi[mix.speciesIndex("N")] = 0.0;
rhoi[mix.speciesIndex("N2")] = 0.78 * rho_init;
rhoi[mix.speciesIndex("O")] = 0.0;
rhoi[mix.speciesIndex("O2")] = 0.22 * rho_init;
double T[2] = {T_init, Tv_init};
mix.setState(rhoi, T, 1); // set state using {rho_i, T, Tv}
// 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_vib(mix);
printResults_vib(0.0, mix);
// Integrate the species mass conservation equations forward in time
double dt = 0.0;
double time = 0.0;
double tout = 1e-6;
while (time < 10.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-9 * std::abs(rho_init * mix.Y()[i] / wdot[i]);
dt = 1e-5;
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_vib(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;
}
/// [example_code]which then gives:
Time[s] T[K] Tv[K] P[Pa] rho[kg/m^3] e[J/kg] Y_N Y_O Y_NO Y_N2 Y_O2
0 4000 200 4341.61 0.00376 2.5854e+06 0 0 0 0.78 0.22
1e-06 4000 117.882 4341.61 0.00376 2.5854e+06 1.46378e-51 3.23063e-26 0 0.78 0.22
0.000101 4000 117.882 4341.61 0.00376 2.5854e+06 2.81548e-29 3.22742e-26 5.06507e-27 0.78 0.22
0.000201 4000 117.882 4341.61 0.00376 2.5854e+06 2.81548e-29 3.22742e-26 1.01903e-26 0.78 0.22
0.000301 4000 117.882 4341.61 0.00376 2.5854e+06 2.81548e-29 3.22742e-26 1.53156e-26 0.78 0.22
0.000401 4000 117.882 4341.61 0.00376 2.5854e+06 2.81548e-29 3.22742e-26 2.04409e-26 0.78 0.22
0.000501 4000 117.882 4341.61 0.00376 2.5854e+06 2.81549e-29 3.22742e-26 2.55661e-26 0.78 0.22
0.000601 4000 117.882 4341.61 0.00376 2.5854e+06 2.81549e-29 3.22742e-26 3.06914e-26 0.78 0.22
0.000701 4000 117.882 4341.61 0.00376 2.5854e+06 2.81549e-29 3.22742e-26 3.58167e-26 0.78 0.22
0.000801 4000 117.882 4341.61 0.00376 2.5854e+06 2.81549e-29 3.22742e-26 4.09419e-26 0.78 0.22
0.000901 4000 117.882 4341.61 0.00376 2.5854e+06 2.81549e-29 3.22742e-26 4.60672e-26 0.78 0.22
0.001001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81549e-29 3.22743e-26 5.11925e-26 0.78 0.22
0.002001 4000 117.882 4341.61 0.00376 2.5854e+06 2.8155e-29 3.22743e-26 1.02445e-25 0.78 0.22
0.003001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81551e-29 3.22744e-26 1.53698e-25 0.78 0.22
0.004001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81551e-29 3.22745e-26 2.04951e-25 0.78 0.22
0.005001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81552e-29 3.22746e-26 2.56205e-25 0.78 0.22
0.006001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81553e-29 3.22747e-26 3.07458e-25 0.78 0.22
0.007001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81554e-29 3.22748e-26 3.58712e-25 0.78 0.22
0.008001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81555e-29 3.22749e-26 4.09965e-25 0.78 0.22
0.009001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81555e-29 3.2275e-26 4.61219e-25 0.78 0.22
0.010001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81556e-29 3.22751e-26 5.12473e-25 0.78 0.22
0.020001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81564e-29 3.2276e-26 1.02502e-24 0.78 0.22
0.030001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81572e-29 3.22769e-26 1.53758e-24 0.78 0.22
0.040001 4000 117.882 4341.61 0.00376 2.5854e+06 2.8158e-29 3.22778e-26 2.05016e-24 0.78 0.22
0.050001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81588e-29 3.22787e-26 2.56275e-24 0.78 0.22
0.060001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81596e-29 3.22796e-26 3.07536e-24 0.78 0.22
0.070001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81604e-29 3.22806e-26 3.58798e-24 0.78 0.22
0.080001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81612e-29 3.22815e-26 4.10061e-24 0.78 0.22
0.090001 4000 117.882 4341.61 0.00376 2.5854e+06 2.8162e-29 3.22824e-26 4.61326e-24 0.78 0.22
0.100001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81628e-29 3.22833e-26 5.12593e-24 0.78 0.22
0.200001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81708e-29 3.22924e-26 1.02534e-23 0.78 0.22
0.300001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81787e-29 3.23016e-26 1.53823e-23 0.78 0.22
0.400001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81867e-29 3.23107e-26 2.05126e-23 0.78 0.22
0.500001 4000 117.882 4341.61 0.00376 2.5854e+06 2.81947e-29 3.23199e-26 2.56444e-23 0.78 0.22
0.600001 4000 117.882 4341.61 0.00376 2.5854e+06 2.82027e-29 3.2329e-26 3.07776e-23 0.78 0.22
0.700001 4000 117.882 4341.61 0.00376 2.5854e+06 2.82106e-29 3.23382e-26 3.59123e-23 0.78 0.22
0.800001 4000 117.882 4341.61 0.00376 2.5854e+06 2.82186e-29 3.23473e-26 4.10485e-23 0.78 0.22
0.900001 4000 117.882 4341.61 0.00376 2.5854e+06 2.82266e-29 3.23564e-26 4.61861e-23 0.78 0.22
1 4000 117.882 4341.61 0.00376 2.5854e+06 2.82346e-29 3.23656e-26 5.13251e-23 0.78 0.22
2 4000 117.882 4341.61 0.00376 2.5854e+06 2.83143e-29 3.2457e-26 1.02795e-22 0.78 0.22
3 4000 117.882 4341.61 0.00376 2.5854e+06 2.83941e-29 3.25484e-26 1.54411e-22 0.78 0.22
4 4000 117.882 4341.61 0.00376 2.5854e+06 2.84738e-29 3.26398e-26 2.06172e-22 0.78 0.22
5 4000 117.882 4341.61 0.00376 2.5854e+06 2.85536e-29 3.27312e-26 2.58078e-22 0.78 0.22
6 4000 117.882 4341.61 0.00376 2.5854e+06 2.86333e-29 3.28227e-26 3.10129e-22 0.78 0.22
7 4000 117.882 4341.61 0.00376 2.5854e+06 2.8713e-29 3.29141e-26 3.62325e-22 0.78 0.22
8 4000 117.882 4341.61 0.00376 2.5854e+06 2.87928e-29 3.30055e-26 4.14666e-22 0.78 0.22
9 4000 117.882 4341.61 0.00376 2.5854e+06 2.88725e-29 3.30969e-26 4.67153e-22 0.78 0.22
10 4000 117.882 4341.61 0.00376 2.5854e+06 2.89523e-29 3.31883e-26 5.19784e-22 0.78 0.22
Equilibrium mass fractions at 4000 K and 4341.61 Pa:
N: 0.00411576
O: 0.223216
NO: 0.0131855
N2: 0.756812
O2: 0.00267114
If I use the one temperature model, I get something that looks way different:
Time[s] T[K] Tv[K] P[Pa] rho[kg/m^3] e[J/kg] Y_N Y_O Y_NO Y_N2 Y_O2
0 4000 4000 4341.61 0.00376 3.36919e+06 0 0 0 0.78 0.22
1e-06 3998.79 3998.79 4340.61 0.00376 3.36919e+06 6.29588e-10 8.15058e-05 0 0.78 0.219918
0.000101 3894.81 3894.81 4253.92 0.00376 3.36919e+06 4.69892e-06 0.0069544 0.000463028 0.779779 0.212799
0.000201 3815.99 3815.99 4186.75 0.00376 3.36919e+06 6.91519e-06 0.0120256 0.00144223 0.77932 0.207205
0.000301 3753.31 3753.31 4132.46 0.00376 3.36919e+06 8.1049e-06 0.0159728 0.00260617 0.778775 0.202638
0.000401 3701.71 3701.71 4087.19 0.00376 3.36919e+06 8.82973e-06 0.0191644 0.00382116 0.778207 0.198798
0.000501 3658.12 3658.12 4048.55 0.00376 3.36919e+06 9.32253e-06 0.0218194 0.00502902 0.777643 0.195499
0.000601 3620.56 3620.56 4014.97 0.00376 3.36919e+06 9.68993e-06 0.024077 0.00620371 0.777094 0.192615
0.000701 3587.68 3587.68 3985.35 0.00376 3.36919e+06 9.9847e-06 0.0260301 0.00733399 0.776567 0.190059
0.000801 3558.52 3558.52 3958.92 0.00376 3.36919e+06 1.02346e-05 0.0277438 0.00841579 0.776061 0.187769
0.000901 3532.4 3532.4 3935.11 0.00376 3.36919e+06 1.04549e-05 0.029265 0.00944856 0.775579 0.185697
0.001001 3508.78 3508.78 3913.47 0.00376 3.36919e+06 1.06545e-05 0.0306287 0.0104335 0.775119 0.183808
0.002001 3351.78 3351.78 3767.16 0.00376 3.36919e+06 1.21145e-05 0.0394214 0.0181595 0.771511 0.170896
0.003001 3262.62 3262.62 3682.3 0.00376 3.36919e+06 1.30609e-05 0.0442339 0.0233087 0.769107 0.163338
0.004001 3202.3 3202.3 3624.27 0.00376 3.36919e+06 1.36984e-05 0.0474463 0.0269473 0.767407 0.158185
0.005001 3157.77 3157.77 3581.19 0.00376 3.36919e+06 1.41294e-05 0.0498169 0.0296082 0.766165 0.154396
0.006001 3123.09 3123.09 3547.53 0.00376 3.36919e+06 1.4418e-05 0.051676 0.0315937 0.765238 0.151478
0.007001 3095.1 3095.1 3520.34 0.00376 3.36919e+06 1.46061e-05 0.0531947 0.0330919 0.764538 0.149161
0.008001 3071.92 3071.92 3497.8 0.00376 3.36919e+06 1.47217e-05 0.054472 0.0342279 0.764008 0.147278
0.009001 3052.32 3052.32 3478.76 0.00376 3.36919e+06 1.47843e-05 0.0555699 0.0350887 0.763606 0.145721
0.010001 3035.5 3035.5 3462.44 0.00376 3.36919e+06 1.48075e-05 0.0565296 0.0357372 0.763303 0.144415
0.020001 2942.09 2942.09 3372.51 0.00376 3.36919e+06 1.42051e-05 0.0623145 0.0369594 0.762733 0.137979
0.030001 2901.25 2901.25 3333.88 0.00376 3.36919e+06 1.34039e-05 0.0652025 0.0356435 0.763348 0.135792
0.040001 2878.27 2878.27 3312.36 0.00376 3.36919e+06 1.2784e-05 0.0669424 0.0343089 0.763972 0.134764
0.050001 2863.93 2863.93 3299.01 0.00376 3.36919e+06 1.23398e-05 0.0680671 0.0332752 0.764455 0.13419
0.060001 2854.52 2854.52 3290.27 0.00376 3.36919e+06 1.2027e-05 0.0688195 0.0325232 0.764806 0.133839
0.070001 2848.18 2848.18 3284.38 0.00376 3.36919e+06 1.18077e-05 0.0693324 0.0319867 0.765057 0.133612
0.080001 2843.83 2843.83 3280.36 0.00376 3.36919e+06 1.16538e-05 0.0696862 0.0316066 0.765235 0.133461
0.090001 2840.83 2840.83 3277.57 0.00376 3.36919e+06 1.15457e-05 0.069932 0.031338 0.76536 0.133358
0.100001 2838.74 2838.74 3275.63 0.00376 3.36919e+06 1.14697e-05 0.0701036 0.0311482 0.765449 0.133288
0.200001 2833.96 2833.96 3271.21 0.00376 3.36919e+06 1.12937e-05 0.070497 0.0307068 0.765655 0.13313
0.300001 2833.81 2833.81 3271.08 0.00376 3.36919e+06 1.12883e-05 0.0705089 0.0306932 0.765661 0.133125
0.400001 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
0.500001 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
0.600001 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
0.700001 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
0.800001 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
0.900001 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
1 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
2 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
3 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
4 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
5 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
6 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
7 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
8 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
9 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
10 2833.81 2833.81 3271.07 0.00376 3.36919e+06 1.12881e-05 0.0705093 0.0306928 0.765661 0.133125
Equilibrium mass fractions at 2833.81 K and 3271.07 Pa:
N: 1.11932e-05
O: 0.0731015
NO: 0.0315333
N2: 0.752352
O2: 0.143002
It would be nice if you could provide an example for a non-equilibrium calculation with the two temperature model. It might only be me, but I think the two-temperature model is what really distinguishes Mutation++ from other packages like Cantera (which can do chemical non-eq.). In particular, with the existing examples I cannot infer how I should update the state for the two temperature model (the update of the state with density and total energy does not seem to give physical results).
In particular, I consider an example along the lines of https://github.com/mutationpp/Mutationpp/blob/master/examples/c%2B%2B/O2_dissociation/O2_dissociation.cpp
which then gives:
If I use the one temperature model, I get something that looks way different: