hy2Foam: Updated turbulent terms and Mach number calculation

applications/solvers/compressible/hy2Foam/createReactingFields.H

@@ -73,3 +73,15 @@ if
  solveSpeciesEqns = true;
 } 
 
+//---------------------------------------------------------
+// Turbulent Schmidt Number - Alex Shepherd
+//---------------------------------------------------------
+
+scalar Sct_
+(
+ readScalar
+ (
+ thermo.transportDictionary().subDict("transportModels")
+ .subDict("diffusionModelParameters").lookup("TurbulentSchmidtNumber")
+ )
+);

applications/solvers/compressible/hy2Foam/createThermoFields.H

@@ -68,10 +68,29 @@ volScalarField Mach
  mesh,
  IOobject::NO_READ,
  IOobject::AUTO_WRITE
+ ), 
+ mag(U)/sqrt(thermo.CpMix()/thermo.CvMix()/psi)
+);
+
+// Changed to include total heat capacities, not just translational-rotational. - Alex Shepherd
+
+volScalarField gamma 
+(
+ IOobject
+ (
+ "gamma",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
  ),
- mag(U)/sqrt(thermo.Cp_t()/thermo.Cv_t()/psi)
+ thermo.CpMix()/thermo.CvMix()
 );
 
+// Write out heat capacity ratio. Usefull for post-processing.
+// Both Mach and gamma fields are updated in upwindInterpolation.H,
+// this wouldn't work otherwise. - Alex Shepherd
+
 volVectorField rhoU
 (
  IOobject

applications/solvers/compressible/hy2Foam/eqns/YEqn.H

@@ -55,7 +55,8 @@ if (solveSpeciesEqns)
  YiEqn += 
  fvm::laplacian
  (
- -speciesDiffusion().rhoD(speciei),
+ - speciesDiffusion().rhoD(speciei)
+ - turbulence->mut()/Sct_,
  Yi,
  "laplacian(rhoD,Yi)"
  )

applications/solvers/compressible/hy2Foam/eqns/eEqnViscous.H

@@ -5,7 +5,11 @@ if (!inviscid)
  fvScalarMatrix eEqnViscous
  (
  fvm::ddt(rho, e) - fvc::ddt(rho, e)
- - fvm::laplacian(thermo.kappa()/thermo.CvMix(), e)
+ - fvm::laplacian
+ (
+ thermo.kappa() / thermo.CvMix() + turbulence->alphat(), 
+ e
+ )
  ==
  fvOptions(rho, e)
  );

applications/solvers/compressible/hy2Foam/numerics/upwindInterpolation.H

@@ -112,18 +112,25 @@ else
  }
 }
 
-//- Fields for boundary conditions
-//volScalarField cp("cp", thermo.CpMix());
-//volScalarField gamma("gamma", cp/thermo.CvMix());
+// - Fields for boundary conditions
+
+// Switched back to using total specific heats for gamma, instead
+// of just translational-rotational. Reference doesn't explain why
+// the choice was made to use gammatr. - Alex Shepherd
+
+volScalarField cp("cp", thermo.CpMix());
+gamma = cp/thermo.CvMix();
+volScalarField c(sqrt(gamma*rPsi));
 
 //- In: G. V. Candler and I. Nompelis.
 // "Computational Fluid Dynamics for Atmospheric Entry"
 // RTO-AVT-VKI Lecture Series 2009-AVT-162, Eq. 3.35, p. 39, 2009
 // + gamma = frozen translational-rotational specific heats of the gas mixture
 // + convention in hy2Foam: translational-rotational specific heats of the gas
 // mixture
-volScalarField gammatr("gammatr", thermo.Cp_t()/thermo.Cv_t());
-volScalarField c(sqrt(gammatr*rPsi));
+// volScalarField gammatr("gammatr", thermo.Cp_t()/thermo.Cv_t());
+// volScalarField c(sqrt(gammatr*rPsi));
+
 Mach = mag(U)/c;
 
 surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);

run/hyStrath/hy2Foam/genericCase/constant/chemDicts/hTCReactionsES

@@ -43,93 +43,101 @@ reactions
  // Reaction no 1
  /*oxygenReaction
  {
- type irreversibleArrheniusReaction; // nonEquilibriumReversibleArrheniusReaction;
- reaction "O2 + O2 = 2O + O2";
+ type nonEquilibriumReversiblethirdBodyArrheniusReaction;
+ reaction "O2 + M = 2O + M";
  forward
  {
- A 7.2e18;
+ A 7.2e15;
  beta -1.0;
  Ta 59340.0;
+ defaultEfficiency 1.0;
  }
  reverse
  {
- A 4.0e17;
+ A 4.0e11;
  beta -1.0;
- Ta 0.0; 
+ Ta 0.0; 
+ defaultEfficiency 1.0; 
  }
  }
 
  // Reaction no 2
  hydrogenReaction
  {
- type irreversibleArrheniusReaction; // nonEquilibriumReversibleArrheniusReaction;
- reaction "H2 + H2 = 2H + H2";
+ type nonEquilibriumReversiblethirdBodyArrheniusReaction;
+ reaction "H2 + M = 2H + M";
  forward
  {
- A 5.5e18;
+ A 5.5e15;
  beta -1.0;
- Ta 51987.0;
+ Ta 51987.0; 
+ defaultEfficiency 1.0; 
  }
  reverse
  {
- A 1.8e18;
+ A 1.8e12;
  beta -1.0;
- Ta 0.0; 
+ Ta 0.0; 
+ defaultEfficiency 1.0; 
  }
  }
 
  // Reaction no 3
  waterReaction
  {
- type irreversibleArrheniusReaction; // nonEquilibriumReversibleArrheniusReaction;
- reaction "H2O + H2O = OH + H + H2O";
+ type nonEquilibriumReversiblethirdBodyArrheniusReaction;
+ reaction "H2O + M = OH + H + M";
  forward
  {
- A 5.2e21;
+ A 5.2e18;
  beta -1.5;
- Ta 59386.0;
+ Ta 59386.0; 
+ defaultEfficiency 1.0; 
  }
  reverse
  {
- A 4.4e20;
+ A 4.4e14;
  beta -1.5;
- Ta 0.0; 
+ Ta 0.0; 
+ defaultEfficiency 1.0; 
  }
  }
 
  // Reaction no 4
  hydroxylReaction
  {
- type irreversibleArrheniusReaction; // nonEquilibriumReversibleArrheniusReaction;
- reaction "OH + OH = O + H + OH";
+ type nonEquilibriumReversiblethirdBodyArrheniusReaction;
+ reaction "OH + M = O + H + M";
  forward
  {
- A 8.5e18;
+ A 8.5e15;
  beta -1.0;
  Ta 50830.0;
+ defaultEfficiency 1.0; 
  }
  reverse
  {
- A 7.1e18;
+ A 7.1e12;
  beta -1.0;
- Ta 0.0; 
+ Ta 0.0;
+ defaultEfficiency 1.0; 
  }
  }
 
  // Reaction no 5
  oxygenAtomicHydrogenReaction
  {
- type irreversibleArrheniusReaction; // nonEquilibriumReversibleArrheniusReaction;
+ type nonEquilibriumReversibleArrheniusReaction;
  reaction "O2 + H = OH + O";
  forward
  {
- A 2.2e14;
+ A 2.2e11;
  beta 0.0;
  Ta 8455.0;
  }
  reverse
  {
- A 1.5e13;
+ A 1.5e10;
  beta 0.0;
  Ta 0.0; 
  }
@@ -142,13 +150,13 @@ reactions
  reaction "H2 + O = OH + H";
  forward
  {
- A 7.5e13;
+ A 7.5e10;
  beta 0.0;
  Ta 5586.0;
  }
  reverse
  {
- A 3.0e13;
+ A 3.0e10;
  beta 0.0;
  Ta 4429.0; 
  }
@@ -157,17 +165,17 @@ reactions
  // Reaction no 7
  waterAtomicOxygenReaction
  {
- type irreversibleArrheniusReaction; // nonEquilibriumReversibleArrheniusReaction;
+ type nonEquilibriumReversibleArrheniusReaction;
  reaction "H2O + O = OH + OH";
  forward
  {
- A 5.8e13;
+ A 5.8e10;
  beta 0.0;
  Ta 9059.0;
  }
  reverse
  {
- A 5.3e12;
+ A 5.3e9;
  beta 0.0;
  Ta 503.0; 
  }
@@ -176,17 +184,17 @@ reactions
  // Reaction no 8
  waterAtomicHydrogenReaction
  {
- type irreversibleArrheniusReaction; // nonEquilibriumReversibleArrheniusReaction;
+ type nonEquilibriumReversibleArrheniusReaction;
  reaction "H2O + H = OH + H2";
  forward
  {
- A 8.4e13;
+ A 8.4e10;
  beta 0.0;
  Ta 10116.0;
  }
  reverse
  {
- A 2.0e13;
+ A 2.0e10;
  beta 0.0;
  Ta 2600.0; 
  }

run/hyStrath/hy2Foam/genericCase/constant/thermoDEM_TRVE

@@ -15,6 +15,16 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+/* Notes
+ * - H2O has 3 rotational dof. Coeff for rot should be 1->1.5
+ * - H2O2 has 3 rotational dof. Coeff for rot should be 1->1.5
+ * - H has incorrect enthlapy of formation ratio. Corrected per NIST
+ * - HO2 has 3 rotational dof. Coeff for rot should be 1->1.5
+ * - NO2 has 3 rotational dof. Coeff for rot should be 1->1.5
+ * - CO vibration mode changed from 2->1
+ */
+
+
 // Atoms and molecules data
 
 N2
@@ -841,7 +851,7 @@ CO
  {
  decoupledCvCoeffs ( 1.5 1 1 0 0 -13290 0 );
  vibrationalList ( 
- 2 3122 
+ 1 3122 
  );
  electronicList ( 1 0 ); 
  }
@@ -917,7 +927,7 @@ H
  }
  thermodynamics
  {
- decoupledCvCoeffs ( 1.5 0 0 0 0 25474 0 );
+ decoupledCvCoeffs ( 1.5 0 0 0 0 26219 0 );
  vibrationalList ( 0 0 );
  electronicList ( 1 0 ); 
  }
@@ -1005,7 +1015,7 @@ H2O
  }
  thermodynamics
  {
- decoupledCvCoeffs ( 1.5 1 1 0 0 -29082 0 );
+ decoupledCvCoeffs ( 1.5 1.5 1 0 0 -29082 0 );
  vibrationalList ( 
  1 2294
  1 5180
@@ -1097,7 +1107,7 @@ HO2
  }
  thermodynamics
  {
- decoupledCvCoeffs ( 1.5 1 1 0 0 1263 0 );
+ decoupledCvCoeffs ( 1.5 1.5 1 0 0 1263 0 );
  vibrationalList ( 
  1 1577
  1 2059 
@@ -1130,7 +1140,7 @@ NO2
  }
  thermodynamics
  {
- decoupledCvCoeffs ( 1.5 1 1 0 0 3993 0 );
+ decoupledCvCoeffs ( 1.5 1.5 1 0 0 3993 0 );
  vibrationalList ( 
  1 930
  1 1900 
@@ -1177,7 +1187,7 @@ H2O2
  }
  thermodynamics
  {
- decoupledCvCoeffs ( 1.5 1 1 0 0 -16393 0 );
+ decoupledCvCoeffs ( 1.5 1.5 1 0 0 -16393 0 );
  vibrationalList ( 
  1 1250
  1 1970